mmam* 


Digitized  by  the  Internet  Archive 
in  2019  with  funding  from 
University  of  Alberta  Libraries 


https://archive.org/details/Sahay1980 


THE  UNIVERSITY  OF  ALBERTA 


RELEASE  FORM 

NAME' OF  AUTHOR  Pratap  Narayan  Sahay 

TITLE  OF  THESIS  SEISMIC  SOURCE  THEORY 

DEGREE  FDR  WHICH  THESIS  WAS  PRESENTED  MASTER  OF  SCIENCE 

YEAR  THIS  DEGREE  GRANTED  Fall  1980 

Permission  is  hereby  granted  to  THE  UNIVERSITY  OF 
ALBERTA  LIBRARY  to  reproduce  single  copies  of  this 
thesis  and  to  lend  or  sell  such  copies  for  private, 
scholarly  or  scientific  research  purposes  only. 

The  author  reserves  other  publication  rights,  and 
neither  the  thesis  nor  extensive  extracts  from  it  may 
be  printed  or  otherwise  reproduced  without  the  author' 
written  permission. 


THE  UNIVERSITY  OF  ALBERTA 


SEISMIC  SOURCE  THEORY 


by 

Pratap  Narayan  Sahay 


A  THESIS 

SUBMITTED  TO  THE  FACULTY  OF  GRADUATE  STUDIES  AND  RESEARCH 
IN  PARTIAL  FULFILMENT  OF  THE  REQUIREMENTS  FOR  THE  DEGREE 

OF  MASTER  OF  SCIENCE 
IN 

GEOPHYSICS 

DEPARTMENT  OF  PHYSICS 


EDMONTON,  ALBERTA 


Fall  1980 


THE  UNIVERSITY  OF  ALBERTA 
FACULTY  OF  GRADUATE  STUDIES  AND  RESEARCH 

The  undersigned  certify  that  they  have  read,  and 
recommend  to  the  Faculty  of  Graduate  Studies  and  Research, 
for  acceptance,  a  thesis  entitled  SEISMIC  SOURCE  THEORY 
submitted  by  Pratap  Narayan  Sahay  in  partial  fulfilment  of 
the  requirements  for  the  degree  of  MASTER  OF  SCIENCE  in 
GEOPHYSICS. 


ABSTRACT 


Theoretical  modeling  of  the  source  region  is  done  in 
order  to  understand  the  physical  process  occurring  in  the 
earthquake.  In  the  beginning  it  was  modeled  as  a  point 
source,  later  on  the  fault  dimension  and  the  slip  on  the 
fault  was  taken  into  account  by  developing  the  dislocation 
models.  Recently  studies  on  the  occurrence  of  an  earthquake 
and  the  resultant  seismic  motion  from  nucleation  of  a  crack, 
its  spreading  and  stopping  which  develops  under  stress 
condition  and  material  properties  of  the  source  region,  have 


been  started. 
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1.  INTRODUCTION 


The  desire  to  achieve  some  depth  of  understanding  of 
the  earthquake  process,  leads  to  the  formulation  of  differ¬ 
ent  models  for  this  complex  phenomenon.  In  seismic  source 
theory,  this  theoretical  modeling  is  studied  systematically. 

The  theoretical  study  of  the  seismic  source  has 
traditionally  been  concerned  with  attempts  to  construct 
models  that  simulate  the  faulting  process,  which,  dating 
from  the  early  work  of  Reid  on  the  great  San- F r ans i sco 
earthquake  of  1906,  is  accepted  as  the  most  valid  physical 
description  of  shallow  earthquakes.  The  first  approach  to 
the  description  of  the  faulting  process  of  an  earthquake  was 
based  on  the  consideration  of  various  simple  point  sources 
and  couples,  within  or  on  the  surface  of  an  elastic  medium. 
Reviews  and  extensive  bibliographies  are  given  by  Balakina 
et  al.  (1961)  and  Honda  (1962)  where  the  original  work  of 
Byerly,  Hodgson,  Nakano,  and  others  is  discussed. 

A  more  elaborate  model  based  on  dislocation  theory  was 
first  proposed  by  Volterra  (1907)  and  later  generalized  and 
extensively  applied  by  Steketee  (1958),  Maruyama  (1963), 
Burridge  and  Knopoff  (1964),  Haskell  (1964,  66,  69),  and 
others . 
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The  modeling  of  the  actual  physical  process  of  the 
source  started  with  Kostrov's  work  on  self-similar  cracks  of 
the  1960's  where  it  is  described  as  crack  propagation. 

With  the  work  of  Burridge  and  Knopoff  (1967)  on  the 
numerical  simulation  of  the  fault  motion  and  earthquake 
occurrence,  computer  modeling  of  the  source  has  started. 

In  chapter  II  and  III  a  closer  look  at  the  theoretical 
nature  of  the  seismic  sources  is  taken.  Chapter  IV  is 
devoted  to  the  earthquake  as  a  fracture  process.  The 
numerical  simulation  of  earthquakes  is  discussed  in  chapter 
V.  Reviews  of  the  source  theory  can  be  found  in 
Archambeau ( 1 968 , 75 ) ,  Ben-Menahem  and  Singh  (1972).  A  summary 
of  later  work  is  given  by  Johnson  (1979). 
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2.  KINEMATICS  OF  THE  SOURCES 


An  earthquake  may  involve  non-linear  processes  such  as 
brittle  fracture,  plastic  flow,  phase  change  and  other 
complicated  phenomena.  As  direct  observation  in  the  region 
where  these  processes  take  place  (source  region)  is 
difficult,  the  problem  is  attacked  by  assuming  some  model 
for  source  region  and  then  calculating  the  displacement 
field  for  such  a  source.  The  objective  of  this  approach  is 
to  estimate  the  parameters  character izi ng  the  source  by  the 
inspection  of  the  radiation  field. 

Inside  the  source  region  the  linear  equations  of  motion 
do  not  hold,  so  the  displacement  field  can  be  calculated  if 
source  region  is  squeezed  to  a  singularity  outside  which 
linear  equations  of  motion  are  valid.  In  such  a  case  the 
equations  of  motion  can  be  solved  by  adding  a  term,  usually 
called,  the  body  force  equivalent  of  the  source.  If  the 
medium  under  consideration  is  isotropic  and  homogeneous  or 
if  the  anisotropy  or  inhomogeneity  is  of  simple  type  then 
one  may  be  able  to  solve  the  equations  of  motion  by  first 
developing  a  Green's  function  with  appropriate  boundary 
conditions  and  then  convolving  it  with  the  source  term. 

Let  ^  be  the  source  term  for  the  earthquake  process, 
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gp 
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then  the  equation  of  motion  are, 


(Cijf<y  -  f  IL  —  '  fc  2.1 

where  CL-^  =  A  kcj  ^  +'^y^|D)  is  the  elastic 

constant  of  the  medium. 

For  homogeneous  isotropic  whole  space,  the  retarded 
Green's  function  for  the  above  equation  with  causality 
cond i t i on  is 

-Yj-Sy)*: 

\ 

.  -1  — I  z  ■ 

.^t-t-H./l5p)-CViYj  -%)^)  S(t-t'-WUO 

where 

X  =  co-ordinate  of  observation  point 

' 

)<  =  co-ordinate  of  source  point 

n.  =  Ix-m 

V L  =  direction  cosine  of  line  from  source  to  obs .  point 

_  CX;  -  *£ 

ft 

Up  =  P  wave  velocity 
\Jb  -  S  wave  velocity 

l  is  the  component  of  displacement  at  the  obs.  point. 

^  is  the  direction  of  the  force  acting  at  the  source. 
Using  symmetry  and  time  translation  property  of  the 


r 

Un 


(f,t;x,t)  =  ^nf)  [(3X 


+Y1yJ(arvj 


' 


If  W': 


5 


Green's  function  we  get 


u.  &  t)  =  Ue  r  g.^5.  ^  y)  i-  *  <o 

J  L  V  5 


•C^K^<V  LljC^  ,  t  )  C^cy 


^  ^  J  JX.  y  t 


2.3 


Putting  2.2  into  2.3  and  carrying  out  the  integration  with 
respect  to  t  ,  we  get 


UL  (2-t )  =  \Qc,  [*,]  A*  *'  +  (  CA<  ^  (kj  [U-hV]  \>k 


\J 


5 


2.4 


i,  ->/ 


■A'tV  4 
5 


where  Guj  is  an  operator  on  <p>  defined  as, 


No: 


y,  -Sy)€3^(.x,t-?,)|  +  yLYj(vviL) 


Vt£ 


^C*'.  t-  VcjO  -%)  (l£x)  '^(x',  t- Vis) 


2.5 


The  above  equation  gives  the  elastic  displacement  at 
any  arbitary  point  inside  a  homogeneous  isotropic  body  V 
enclosed  by  S  in  terms  of  distribution  of  ^  in  V  and 
the  value  of  llj,  and  U ^,cy  on  the  SD  .  Now  splitting  5 
into  an  inner  S0  (surrounding  the  source)  and  an  outer 


6 


boundary  St  of  V  and  letting  5,  to  go  to  infinity  where 
homogeneous  boundary  conditions  are  satisfied,  the  surface 
integration  in  the  equation  2.4  reduces  only  on  S0 


The  various  earthquake  source  models  are  basically 
different  representat i ons  of  the  source  term  ,  the 
source  region  boundary  ^  .  inside  which  non-linear 
phenomena  are  taking  place,  and  the  displacement  and  stress 
fields  on  ^  .All  theoretical  source  representations  can 
be  classified  as : 


1.  Point  source  models 

2.  Dislocation  models 

2.1  Point  source  models:- 


In  the  point  source  representation  of  the  earthquakes 
the  source  region  is  squeezed  to  a  point,  so  that  the 
surface  term  is  expressed  as  a  single  force  or  a  cluster  of 
forces  localized  in  a  squeezed  source  region.  Therefore  the 
expression  for  the  displacement  field  in  this  case  becomes 


U;(x,t)  = 


2.6 


Using  the  above  formula,  the  displacement  field  for  any 


' 


point  source  can  be  calculated. 


Single  Force : - 


The  trivial  representation  for  a  point  source  is  a 
single  force.  The  source  term  for  a  single  force  is 


f.  =  e-  Sy  •  F(t) 


2.7 


fit)  is  the  time  dependent  part  of  the  force.  With  this  we 
get , 


2.8 


-  C ^  '  £y 


0- 

6' 


FCt-  *-/ir.) 


The  first  term  in  the  equation  2.8  attenuates  like  Tc 
and  is  called  the  near  field  term  while  other  two  terms 
which  attenuate  like  K  are  called  the  far  field  terms. 


Near  field  term:- 


8 

This  can  be  decomposed  into  the  terms  apparently 
propagating  with  P  and  S  waves  velocities  by  introducing  the 
following  integral, 

t 

Tct)  = 

o 

Integrating  by  parts, 

(3  FCt-3)  =  -  I  1  Ct  -  +  &  i(t  -Vtrf ) 

Ki$? 

-  T  (t-Vir*)  +  T(t-^ArP). 


FCn.H'l 


2.10 


Putting  2.11  into  the  equation  2.9,  we  find  four  terms  like 

j^TCt-VujO,  ^TCt-V*),  2. 12 

The  first  term  propagates  with  P  wave  velocity  and 
attenuates  as  II3  ,  with  the  wave  proportional  to  a  double 
integration  of  the  force.  Likewise,  the  third  term  also 
propagates  with  P  wave  velocity  and  attenuates  as  ,  with 
the  wave  proportional  to  one  integration  of  the  force. 
Similarly  other  two  terms  propagate  with  S  wave  velocity. 
This  theoretical  decomposition  is  elegant,  but  on  a  record 
all  arrive  at  the  same  time  or  nearly  at  the  same  time  so 
they  cannot  be  separated.  They  can  be  distinguished  only  if 
they  are  observed  at  a  long  distance.  But  due  to  the  fast 
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attenuation  with  distance  they  are  negligible  there. 
The  radial  component  of  u!J  is, 

Mr* 

'(ill"  = 


Therefore  the  S  wave  in  this  form  will  have  a  radial 
component . 

Define  a  vector  Y-7  orthonormal  to  Y* 


N 

Then  the  transverse  component  of  Lie  is, 

n/\y6 

Hi-  U 1  =  -  Wtt  f)  Vj 

V(f 

Clearly  the  P  wave  here  will  have  a  transverse 
component . 

F ar  field  term : - 

The  far  field  term  has  two  parts,  one  propagates  with  P 
wave  velocity  while  another  with  S  wave  velocity,  both  are 
proportional  to  the  applied  force  and  attenute  as  .  The 

far  field  P  term  is, 


l  F(t-1)d| 


2.15 


. 
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/  ? 

As  V L'U-i~0  ’  the  ^ar  field  P  term  does  not  have  a 

transverse  component.  The  radial  component  is 


y.  uj  =  (4irf)'  Yj  (^)  F(t->VirP)  2.16 

Clearly  the  directivity  pattern  is  proportional  to  ^  , 
the  direction  cosine  of  line  joining  source  and  observation 
points,  looks  like  fig.  2.1a  . 


*  -  0Wcy 

-  =  toward 


a 

P  wave 

Fig.  2 . 


b 

S  wave 

Directivity  patterns  for  single  force 


The  far  field  S  term  is 


=  - ^ n  f )  (XCV,  FCt-V^ 


2.  17 


The  radial  component  is  zero  while  the  transverse 
component  is 

tf-U-  =  (Snf)  V  FCt  -V*) 


2.  18 


Here  the  radiation  pattern  is  proportional  to  Y:  ,  and 

r 

looks  like  fig.  2.1b  . 

Single  coup  1 e : - 

The  source  term  for  a  single  couple  due  to  forces 
• 

acting  in  the  y  direction  separated  along  the 
d i rect ion  is, 


If  the  source  term  is  called  a  single  couple 
without  moment  or  dipole,  otherwise  single  couple  with 
moment  (fig.  2.2).  Thus  there  can  be  three  dipoles  and  six 
coup  1 es . 


Fig.  2.2a  Dipole  cores spound i ng  to  (jk)=(11),  (22),  (33) 


Fig.  2.2Ds,ngie  couple  (23)  and  (32). 


rfcrt 
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Now  the  source  term  for  a  single  force  (equation  2.7) 
on  differentiation  w.r.t.  yields  the  source  term  for  a 

single  couple  (equation  2.19),  therefore,  the  solution  for  a 
single  force  (equation  2.8)  on  differentiation  w.r.t.  XK 
will  give  the  solution  for  a  single  couple,  which  is 
(Haskell  1963), 


(tmf)'  [(sC^k  +  ck  is 


+  -TCt-Woi)/^4 


2.20 


The  directivity  pattern  in 
and  a  single  couple  with  moment 


the  far  field  for  a  dipole 
i s  shown  in  fig  2.3  . 


S  wave 


Fig.  2.3a  Directivity  patterns  for  doublet. 
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* 


toword 


P  wave 

5  wave 

Ficj.  2.3b  Directivity  patterns  for  single  couple  witt 
moment  s 


Center  of  dl 1  a  tat  ion ,  Doub 1 e  coup  1 e  and  Center  of  rota  t i on : - 

The  source  terms  for  these  are  linear  combinations  of 
single  couple  source  terms. 

The  center  of  dilatation  is  a  linear  combination  of 
three  equal  and  mutually  perpend i cu 1 ar  dipoles. 

The  double  couple  is  linear  combination  of  two  coplanar 
single  couples  of  equal  but  opposite  moments  with  mutually 
perpend i cu 1 ar  force  directions. 

The  center  of  rotation  is  a  combination  of  two  coplanar 
single  couples  of  moment  of  equal  magnitude  acting  in  the 
same  direction,  so  that  the  net  moment  is  not  zero  but  the 
net  force  is  zero.  Although,  this  type  of  source  is  not 
possible  in  the  real  earth. 

The  distribution  of  forces  for  all  of  these  are  shown 

in  fig.  2.4  . 


(12)  +  (21) 


Fig.  2.4a  Three  fundamental  double  couples 


Fig  2.4b  Center  of  rotation  with  moment  about  X(  ,  ,  and 

-d i rect i ons . 


Clearly  the  displacement  field  for  each  can  be  obtained 
with  a  proper  combination  of  solution  for  single  couples 
( equa t ion  2.20). 


The  directivity  pattern  for  each  is  shown  in  fig.  2.5 
3 

1 


Fig.  2.5b  Directivity  patterns  for  double  couple 


' 
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As 

var ious 

source 

terms  are 

just 

combinations  of 

single 

couple 

terms , 

thi  s 

leads  to 

the 

ambiguity  in 

the 

i nterpretat i on . 

The 

combi nat i on 

of 

two  oppos i te ly 

signed 

dipoles 

at  right 

angles  to  each 

other 

is  the  same 

as  a 

double  couple  (fig.  2.6a).  Similarly  a  suitable  combination 
of  a  double  couple,  a  center  of  dilatation  and  a  linear 
doublet  is  equivalent  to  a  differently  oriented  double 
couple  (fig.  2.6b). 


(O 


F > g .  2  6  Ambiguity  in  the  source  representation. 
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Table  2.1  The  source  term  for  var i ous  sources 
Source 

f- 


Single  force 

1  . 

in  the 

x^-d i  rect  i  on 

A 

£ 

2. 

i n  the 

x-di rect i on 

A 

7- 

3. 

in  the 

x-di rect i on 

A 

e: 

Dipole 

A 

1  . 

in  the 

x-di rect i on 

-e, 

2. 

in  the 

x^-di  rect  ion 

\ 

-  e. 

3. 

in  the 

x3-di  rect  ion 

A 

-£. 

Single  couple 

1  . 

(  12) 

A 

-e, 

2. 

(21  ) 

A 

-<7 

3. 

(23) 

A 

"^2. 

4  . 

(32) 

A 

-e. 

5. 

(31  ) 

A 

-  e: 

6. 

(13) 

A 

-e, 

-e,  Mx, -a,')  Scxi-xao  ^<7*3 

^x 

5  C.X,-X/ )^&  (-*-  *  _:t2j 
x  xxo. 


x4) 


3  SCx-X,')  *3 

^  (33-3 


Xi. 


iPr^1)  s c^p'J 


3 


^  S(p-x;)  ^  o  3; 

>3-1 


I  fcOk-IV 


h 


Doub 1 e  coup  1 e 

23 ) ♦  ( 32 1  - [X  ec*.-*^ s ^3 asc^^O sci 3-^ 


1  . 
2. 
3. 


~  £"  7"  777  a  ')  X  §77x  a 

(  3 1 )  +  ( 1 3 )  -  )  OLX-I-X})  X 

',  ,|,-n  C,  -,/)+£  3  §OlrP|',)  SCls-^dX 

1 1 2 1  +  ( 2 1 )  -te, 


Center  of  rotation 

1.  (321-123)  'i^fuUl'^n^XX'CtU^^^X\  ,n,l  X’} 

h  _e,  S<tx  i  - ~x(J&  C 3] bCxy-fi) 

2sX,  "  ‘  '  2X2 


2.  (  131-131  )  -IX  &w'  •  dU 

3.  ( 2  1  )  -  (  1 2  I  -L'^.?i7W'';t',')^W*’’;I,J 

Center  of  compress i on 


- 


)  &(Xv-x7>  SCXi'^72  ^a.  6Ur*,)|- ' 


A 


3!  )  ^  ^3 


6t^-X,7  feUa-ia  BSUj-Xjj] 

"5  X-j 


. 
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The  source  term  fc  for  various  point  sources  are 
presented  in  the  table  2.1. 

The  accepted  model  for  the  shallow  focus  earthquake  is 
a  double  couple.  Haskell  (1963)  has  given  directivity 
pattern  for  a  double  couple  in  half  space  for  Rayleigh  wave. 
Burridge  et  al  (1964)  and  Alterman  et  al  (1970)  have  studied 
the  effect  of  burying  a  double  couple  directly  beneath  the 
free  surface  of  an  elastic  half  space. 

For  the  deep  focus  earthquakes  Knopoff  and  Randall 
(1970)  proposed  a  combination  of  double  couple  and  center  of 
dilatation  called  compensated  linear  vector  dipole  (fig. 
2.6c  ).  This  combination  does  not  produce  volume  change,  has 
no  net  force,  no  net  moment  and  has  five  degrees  of  freedom 
in  spatial  orientation. 

2.2  Dislocation  model 


In  the  dislocation  representation  of  the  earthquakes 
the  source  region  is  shrunk  such  that  the  source  boundary 
-*0  may  be  regarded  as  madeup  of  two  sides  f>D  and  SD  .The 
displacement  and  stress  fields  on  each  surfaces  are 


Fig.  2.7  Dislocation  surface. 
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( u/  - U-j  )  .  the 


considered  to  be  different.  Therefore  in  the  surface  term  of 
the  equation  2.3,  we  get  a  term  like 
superscript  +  refers  to  the  value  at  Si  and  -at  s~,  because 
the  normals  to  the  surface  are  opposing  each  other. 

Def i ni ng 

C^-u;)=  auj  22i 

which  may  be  regarded  as  a  jump  in  the  displacement  and 
stress  fields  on  S0,  and  putting  it  into  the  eq .  2.3,  we  get 


-kj 

f 


+-0 


U-:(x /U  - 


it' [Gcj  Cx,t;  +  AUf(f,t'jCi 


-*Q  V 


'  S  Cj  (>C-5t  )1/t  )  ?t 


2.22 

a-? 


In  the  above  equation  for  the  surface  integration  a  new 
variable  is  used  to  distinguish  as  it  runs  only  on  the 
surface  5  ,  and  hereafter  S0  will  be  cal  led  ]>  . 


Introducing  the  Radon  transformed  Green's  function 

,3-w 


Gcj  (2  ,t ;  f ,  f  )  =  ^  (£  f  )  &y  (*,  t ,  ?,  t'U  1 


V 


V 


f  ~y 


2.23 


plugging  it  into  the  equation  2.22,  we  get 


+*0 


-  \J  ^  J 

+  &awcf  ,t'j  p/ty1' 


-> '  %  \ 
*-X) 

2.24 


' 


■ 
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The  solution  with  homogeneous  boundary  conditions  on 
5,  without  any  jump  conditions  on  TE.  is  simply 
space-time  convolution  of  the  Green's  function  with  body 
force  term.  Since  in  the  above  equation,  the  surface 
integral  within  middle  bracket  is  involved  in  the  same  way 
as  £  ;  it  may  be  concluded  that  the  effect  of  a  jump 

across  .  has  the  same  effect  as  introducing  an  extra 

body  force  term  ,  given  by 


aO)  =  -  +  ^,vtXb)  ^5 


2.25 


into  an  unfaulted  medium.  Therefore  the  Green's  function  for 
an  infinite  domain  is  valid  here  too. 

Us  i  ng 

aTj,  -  \  ^ 

we  get  a  relation  expressing  the  equivalence  of  the  body 
force  and  the  discontinuity  of  displacement  and  stress  as 

=  -  \  +  2.26 

Putting  2.26  into  2.24  and  performing  the  time  integration 
we  get  the  displacement  field 

uc(.£t)  - 

v/ 


2.27 


* 
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where  is  defined  as  in  the  equation  2.5  . 

As  is  a  body  force  in  the  source  free  region,  so 

neglecting  it  we  get  the  displacement  field  due  to  the 
dislocation  source 

Hi[x,  0  -  -  j 

Thus  Knowing  the  dislocation  surface  or  source  boundary 
^  and  the  jump  discontinuity  across  it,  the  body  force 
term  can  be  ca lcul a  ted,  which ,  on  integration  with  the  full 
space  Green's  function,  will  give  the  displacement  field. 


In  an  earthquake  on  the  dislocation  surface,  the  stress 
discontinuity  does  not  exist  because  along  the  dislocation 
direction  the  component  of  stress  vanishes  and  is  continuous 
in  all  other  directions.  Therefore,  the  term  may  be 

"V  V 

dropped  in  the  equation  2.26  and  AUj  (3  A )  will  be  the 
displacement  field  jump  across  Tel  .  The  expression  for  the 
far  field  P  wave  can  be  obtained  using  the  Green's  function 
for  that  part  only  as 


u[(*,t)=  -  OTS?) 


2.29 


Am,  0*  » 

2.30 


. 
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In  the  equation  2.30  another  three  terms  attenuating 

-2-  -*•  N 

as  fa.  are  dropped.  Further,  in  the  far  field  tz.  =■  \  *  -  *%  \ 
is  much  larger  than  the  dimension  of  Z  ,  so  all  terms 
involving  fa  in  the  equation  2.30  may  be  regarded  as  nearly 

constant  and  can  be  pulled  outside  the  integration,  except 

$ 

the  time  delay  term  in  AU-J,  which  will  vary  strongly  for 
small  changes  of  fa  .  Assuming  the  dislocation  is  linearly 
polarized  in  the  direction  Ylj.  ,  ie,  Atlj, ALL  C3*t)  we 
have 


U ■  (A -to  x  nf  fa)  '  (  Cj.K t>«v' v‘  V Yv ) (■ i u-  C? . tr ■ - tyirf )  If 


2.31 


r 


Similarly  the  far  field  S  wave  term  is: 


It  is  evident  that  the  P  wave  is  only  longitudinal 
while  the  S  wave  is  only  transverse.  The  first  part  in  the 

above  equations  2.31  and  2.32  suggests  that  both  P  and  S 

-l 

waves  attenuate  as  %  .  The  second  part  involving  the 
direction  cosine  of  the  point  of  observation,  the 
dislocation  surface  or ientation( and  the  direction  of  the 
dislocation  )  gives  the  directivity  pattern.  The  third 

part  involving  integration  of  the  discontinuity  over  the 
dislocation  surface  gives  the  wave  form  and  has  information 
about  the  dislocation  surface  dimension  and  the  amount  of 
dislocation.  In  observational  seismology,  an  idea  about  the 


' 
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orientation  of  the  dislocation  surface  and  dislocation 
direction  is  obtained  from  the  directivity  pattern.  By 
analysing  the  wave  form  (actually  spectra),  an  estimate 
about  fault  plane  dimension  and  the  amount  of  dislocation  is 
obtained . 

Let  us  examine  the  directivity  pattern  for  a  slip 
di s locat ion . 

Here  ,  so  we  have  only  two  non  vanishing  terms 


)p=  K  ,  ,  and  ,  ‘V-K 


and  ^ 


For  the  P  wave  the  directivity  pattern  yields  si  X i  Yj  Y# 


and  for  the  S  wave  ju.\  (  YL  iT  -  $L'j  )  4 Yj  (  XK  -  j 

These  terms  are  same  as  the  radiation  pattern  for  far 
field  P  and  S  terms,  respectively,  for  a  double  couple. 
Since  we  do  not  have  a  priori  Knowledge  of  dislocation 
surface,  it  may,  to  a  first  appoximation  be  shrunk  to  a 
point,  ie, 


All  Sis,)  U\)  F(t) 


SO 


which  is  same  as  the  wave  term  for  double  couple.  Thus  in 
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all  respects  a  point  slip  dislocation  is  equivalent  to  a 
double  couple  in  an  unfaulted  medium. 

For  any  arbitrary  dislocation  surface,  one  may  create  a 
source  such  as  a  distribution  of  point  sources  over  a  finite 
region  acting  simultaneously.  Alternatively  one  can  have  a 
point  source  tracing  out  a  moving  pattern  over  a  finite 
volume  of  material.  The  first  case  can  be  thought  of  as  the 
infinite  velocity  limit  of  the  second  case.  In  either 
example,  one  has  to  formulate  the  point  source  problem  from 
the  point  dislocation  (surface)  model,  and  then  sum  up  the 
effect  of  the  distribution  by  integrating  over  space. 

Ben-Menahem  (1961)  studied  the  radiation  pattern  of 
surface  waves  and  later  (Ben-Menahem,  1962)  for  body  waves. 
For  a  simple  model  of  fault  by  considering  a  moving  point 
source,  he  has  shown  that  the  finiteness  of  the  seismic 
source  plays  an  important  role  in  the  radiation  pattern 
whenever  the  dominant  wavelength  of  the  source  spectra  is  of 
the  order  of  the  dimension  of  the  fault  or  when  the  time  of 
rupture  is  of  the  order  of  the  period.  Here  the  spectrum  is 
modulated  by  5tv-x/x  type  of  function, 


A  (JjT,  6 ) 


(9  =Azimuthal  angle  of  point  of  observation  from  source 


=Source  length 


; 


. 
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Vs-  Velocity  of  propagation  of  particular  phase 

=veloci ty  of  propagation  of  moving  source  or  rupture 
velocity 

It  was  also  demonstrated  that  effective  fault  length 
and  rupture  velocity  can  be  recovered  by  taking  the  ratio  of 
the  displacement  of  the  same  spectral  component,  at  two 
stations  which  are  located  symmetrically  with  respect  to  the 
source.  This  ratio  is 


the  vector  D  is  called  vector  directivity;  its  magnitude 
and  phase  are  obtained  for  all  spectral  frequencies  and  by 
comparing  it  with  theoretical  results,  the  fault  parameters 
b  and  v*,  are  deduced. 


Hirasawa  and  Stauder  (1965)  and  Savage  (1965) 
considered  that  the  time  after  which  the  wave  arrives  at  the 
observation  point  is  short  compared  to  the  time  for  rupture 
to  traverse  the  source  region.  Thus  they  ignored  the  effect 
of  cessation  of  motion  (stopping  phase).  Savage  explicitly 
gave  the  following  modulating  function  in  this  case 
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for  unilateral  faults  and 

2V,, 

Fa®)  —  ; - a  7 

1 1  -  (V^/vF)  UnO  j 

for  bi lateral  f au Its. 

This  yields  that  the  nodal  planes  of  the  original 
directivity  pattern  are  unchanged;  a  point  made  by  Hirasawa 
and  Strauder. 

Savage  has  also  shown  that  the  resulting  response  is 
roughly  rectangular  and  has  time  duration  proportional  to 
b /  F,(e)  •  Therefore  the  product  of  amplitude  (which 

depends  upon  F|(0J  and  the  time  duration  does  not  depend 
strongly  upon  the  direction  of  dislocation  propagation  since 
the  factor  Ft(£) )  is  cancelled.  This  is  essentially  same  as 
Ben-Menahem  result  for  small  b  or  large  I  . 

The  analysis  for  arbitrary  source  can  be  systematized 
by  expanding  the  Green's  function  in  vector  spherical 
harmonics.  Carrying  out  the  integration  2.3  leaves  the 
displacement  represented  as  a  summation  over  vector 
spherical  harmonics,  a  multipole  expansion.  One  then 
identifies  the  individual  terms  with  simple  forces,  couples, 
double  couples,  center  of  dilatation,  etc.  (  Archambeau 
(1968),  Ben-Menahem  and  Singh  (1972),  Randall  (1972)).  An 
extensive  study  on  unified  approach  to  the  representation  of 
seismic  sources  has  been  done  by  Singh  et  al  (1973). 


' 
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A  more  sophisticated  approach  is  developed  by  Backus 
and  Mulcay  (1976a, 76b) ,  and  Backus  (1977a, 77b)  who  presented 
a  thorough  analysis  of  the  sources  in  tensor  formulation. 
They  describe  a  number  of  different  ways  of  representing  an 
arbitrary  indigenous  seismic  source,  introduce  the  concept 
of  stress  glut,  discuss  the  uniqueness  of  the  various  source 
descriptions,  and  outline  a  method  of  expanding  a  source 
representation  in  a  series  of  polynomial  moments. 

A  more  realistic  model  for  an  arbitrary  source  will  be 
to  assume  that  the  point  dislocation  propagates  itself 
instead  of  the  point  source  equivalent  of  point  dislocation. 
Here  a  rupture  criterion  for  the  propagation  of  a  point 
dislocation  can  be  incorporated  to  take  into  account  rock 
strength,  friction  etc.  This  approach  is  known  as  crack 
propagation.  In  the  past  few  years  extensive  work  has  been 
done  on  this  aspect  of  dislocation  model.  In  chapter  IV  we 
will  systematically  outline  the  development  of  this 


approach . 


. 


3.  SEISMIC  SPECTRUM 


We  Know  that  a  waveform  is  given  (say  for  P  wave)  by 


3.  1 


Therefore  the  observed  waves,  ie,  seismograms,  can  tell  us 
about  the  fault  dimension  2T  and  dislocation  A  . 

A  practical  approach  to  this  inverse  problem  is  to 
describe  7Z  and  ZU,  with  a  small  number  of  parameters.  We 
consider  the  simplest  of  such  models:  a  unidirectional 
propagating  dislocation. 

We  assume  that  ZL  is  rectangular  with  the  longer  side 

as  fault  length  (  L  )  and  shorter  side  as  fault  width  ( U/  ) . 

— > 

We  also  assume  that  the  dislocation  starts  at  f  ~ O  and 
propagates  along  the  fault  length  with  rupture  velocity  [y  , 
and  that  the  dislocation  has  identical  form  at  any  .  We 
call  this  time  function  7  source  time  function'  and  write  it 

D(t)  • 

By  the  above  parameterization  AU.(3  O  can  be  wpi^en 
as 

AlUf  ,0  -  D(t  -Si  /y) 

3.2 

where  the  J  axis  is  taken  along  the  fault  length  L  . 
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I  '■'?  f1!'  Up  tie: 


. 


28 


With  above,  the  expression  3.1  becomes 

=  uj  ji)(t-H/cf-?,/i 3.3 
r  0 

This  gives  the  wave  form  for  far-field  P  waves.  Since  we  are 
concerned  about  the  far-field,  we  can  put 

::r  Jl-o  ^  CcoQ 


Fig.  3.1 


Then  the  equation  3.3  may  be  written  as 

L 


UJ  L  f  •  . 


o 


3.4 


t  -  ?!  '/lt  -  ^©/irf ) 

rL  =L  L'/v-W%) 


where 


. 


29 


The  above  equation  shows  that  the  observed  wave  form  is 

0 

a  moving  average  of  D  over  the  time  interval  ~T^ 


From  the  equations  2.31  and  3.4  the  far-field  P  waves 
can  be  expressed  as 

I 

F(^)  =Yl-Ulp(x,t)=(AnfirfV)  m-Kfo-t'M&.s 


L 

o 


In  order  to  determine  the  source  parameters  by 
comparing  the  above  theoretical  wave  form  with  the  observed, 
it  is  convenient  to  work  in  the  frequency  domain.  We  define 
'the  seismic  spectrum'  as  the  Fourier  transform  of  the 
ground  displacement: 


f 


a  , 


ur  1 


=  ^  Wy1j 


UjL 

\ 


/  —uUj  t 


h 


/0 


J 

o 


Defining 


we  have  , 


4^3 

f  , 

Dlt)  c  dUt 


3.6 


e  vv) 


-LorVi 

& 


Here  the  amplitude  spectrum  is 


3.7 


•  sq  s:ii  :<  sn.  9.>-m-?rsfc  -o.  qjt  1,3  r; 


and  phase  spectrum 

__  1  yi  4-  TI  -h 

\TP  * 

where  ^(urj  is  the  phase  spectrum  for  dislocation  ve loci  ty  \J(ur) 

,  -  i  ^4) 

VJat)  -  |  vVj  |  e 

The  equation  3.7  describes  the  charecter i s t i cs  of 
amplitude  spectra  of  P  waves.  Similar  expression  may  be 
obtained  for  S  waves.  Let  us  examine  the  spectrum  starting 
from  very  low  frequency.  For  small  \jX  ,  X  is  small  so 
SwiX/V  is  nearly  unity.  Also  |  V(ur)(  ,  for  small  UX"  ,  is 
almost  a  constant,  equal  to  permanent  off-set.  From  equation 
3.6  for  zero  frequency 

V[v)  -  jl)(- t)  di  -  I )(>^)  - D(-°o  ) 

"/0 

which  is  equal  to  permanent  off-set  A  D  across  the  fault 
plane.  Thus  from  the  equation  3.7,  we  have 

r ( >?, 0)  =  (Anfl^X  wLaj 

Since  the  maximum  value  of  V^y^v^flj  is  /t 

for  slip  dislocation  in  an  isotropic  body  then 

'Ffio)  =  id L  ad 


■ 
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The  quantity  on  the  right  hand  side  is  called  seismic 
moment,  can  be  readily  obtained  from  the  zero  frequency 
amplitude.  It  is  also  evident  that  the  spectrum  at  low 
frequencies  will  be  flat. 

For  higher  frequency,  the  factor  SuyiX/X  i  n  the 
equation  3.7  produces  a  B  dependent  spectrum.  The  envelope 
of  $uaK./)(  is  proportional  to  UT  .  This  smoothing  effect 
is  weakest  in  the  direction  of  rupture  propagation  (  9-0) 
and  strongest  in  the  opposite  direction  (  IT).  As  a  result 
more  high  frequency  waves  will  be  observed  in  the 
propagation  direction. 


The  \J(yr)  part  will  also  modulate  the  spectrum 
according  to  the  choice  of  source  time  function  D(tX  .  The 
simplest  form  of  the  function  may  be  a  step  function.  In 
that  case  is  proportional  to  and 

\fo)  -  (  i(t)  it  =  GsmUt 

So  the  step  function  is  physically  impossible  for 
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spontaneous  rupture,  because  the  far-field  wave  will  have 
infinite  seismic  energy.  Therefore  'rise  time'  ,  a  finite 
duration  over  which  1)  (£)  r i ses ( Ben-Menahem  and  ToKsoz 
(1963),  Haskell  (1964)),  is  introduced.  Usually  they  are 
assumed  to  be  ramp  or  exponential  functions. 


[V  \UT)  |  - 


For  a  ramp  function 


For  a  exponential  function 


In  any  case,  the  finite  rise  time  introduces  an 
additional  factor  proportional  to  '  at  high  frequencies. 
Thus,  we  have  a  general  feature  of  the  seismic  spectrum  as 


33 


shown  in  the  fig.  3.3. 


^'9-  3.3  Seismic  spectrum. 


There  are  other  dislocation  models  in  which  the  rupture 
propagation  is  two-dimensional  (  Savage,  1966).  In  such 
cases,  the  smoothing  effect  at  higher  frequencies  is 

-‘L  _  | 

proportional  to  UT  instead  of  ur  of  the  uni -di rect iona 1 
propagating  dislocation.  This  gives  ur3  dependence  at  the 
highest  frequencies. 

Apart  from  the  low  frequencies  constant  level  and 
higher  frequencies  proportionality  to  negative  powers  of 
frequency  there  is  another  feature  in  the  spectrum  which  is 
defined  by  the  intersection  of  the  low  and  high  frequency 
trends.  This  is  called  ' corner  frecquency'  and  is  related  to 
the  length  of  the  fault  (Brune,  1970). 
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For  an  ensemble  of  earthquakes  with  a  wide  range  of 
magnitude  we  can  construct  a  scaling  law  of  seismic 
spectrum.  For  example  ,  i f  we  assume  that  a  large  earthquake 
and  a  small  one  are  physically  similar,  we  find  out  that  the 
source  parameters  should  satisfy  following  conditions, 

^  U  AD  b(T" 


and  \T  =constant . 

Then,  the  seismic  spectrum  will  take  following  shape  for  the 
case  of  UT'  model  (  Aki  ,  1967). 


Period  Iseci 


0002  0.01  005  0.2  05  12  5  10 


Frequency  (Hzl 


Fig.  3.4  Seismic  spectra  for 
an  ensemble  of  earthquakes 
(Aki.  1967  ) 

The  broken  line  is  locus  of 

corner  frequency,  proportional 
-3 

to  UT  . 


. 
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From  the  similarity  assumption,  the  corner  frequency 
inversely  proportional  to  the  source  length  and  seismic 
moment  is  proportional  to  the  cube  of  source  length. 
Therfore,  the  seismic  moment  is  inversely  proprotional  to 
the  cube  of  the  corner  frequency.  The  results  are  shown  in 
the  fig.  3.4. 

Simple  models  discussed  above  are  extensively  used  for 
the  source  parameters  determination.  For  predicting  seismic 
effects  of  dangerous  earthquakes  or  distinguishing  nuclear 
tests  more  complicated  model  should  be  considered. 


4.  DYNAMICS  OF  THE  SOURCES 


In  order  to  find  the  displacement  field  from  a 
dislocation  model  of  the  earthquake  an  intregation  over  the 
dislocation  surface,  TE(x,t)  ,  of  the  displacement  field 
discontinuity,  &Uk(x,t),  has  to  be  carried  out.  As  a  first 
order  of  appoximation  the  inhomogeneities  present  in  the 
earth's  crust  seem  enough  to  cause  variation  in  the  rock 
strength  and  friction  and  to  create  the  local  stress 
build-ups.  A  stress  concentration,  sufficient  to  overcome 
the  resistance  of  the  rock  to  fracture,  will  start  the  crack 
and  when  stress  diminishes  or  the  resistance  increases  the 
crack  will  stop.  Therefore,  it  may  be  visualized  that 
depending  upon  conditions  of  stress,  rock  strength  and 


friction 

a  rupture  front  nucleates 

from  a 

poi nt ,  then 

spreads 

over  a 

region  with  some  finite  velocity  and  finally 

ceases 

to  run 

,  resulting  in  some 

distribution  of  the 

d i sp 1 acement 

field  discontinuity 

over  the 

dislocation 

surface 

formed 

due  to  the  crack. 

Thus  the 

di s locat ion 

surface 

i  s  a 

function  of  space 

and  time 

and  has  to  be 

determined  along  with  the  discontinuity  of  the  displacement 
field  over  it  from  the  condition  of  initial  stresses,  rock 
strength  and  friction.  Mathematically  the  problem  may  be 
formulated  as  follows. 
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Let  the  medium  with  elastic  moduli  and  density 

f  be  held  by  a  constant  stress  Cf  °j  .  At  time  t  =  0  ,  a 
crack  nucleates  at  the  origin  and  spreads,  so  at  time  t>0  it 
covers  the  surface  x(x,  t)  .  Let  Lt((x,t)  be  the 
displacement  field  measured  relative  to  the  state  of  strain 
prevai 1 i ng  for  t=0  , 

We  have, 


A 


hV  ^-1 °rcV 


and  the  total  stress  in  the  medium  is  then 


cf.  • 
LJ 


-  <ruc-  -t-  OX-; 
lJ  lJ 


Due  to  the  crack  there  is  relative  movement  of  the 

A 

parts  of  the  medium  adjacent  to  the  different  sides  of  2_  ; 
if  £  denotes  the  limiting  value  of  Uc  on  two  sides  of 


u, 


-  u,  -  AU, 


If  denotes  the  positive  normal,  then  for  tensile  cracks 


A1A  •  >0 


AU,  X  Y\  -  o 


And  for  sliding  crack 


AUx  •  nc  =  o 


AUU  X  V)',  70 
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The  position  of  area  of  crack  at  each  moment  is 
determined  by  a  contour  surrounding  the  edge  of  the  crack; 
let  be  the  unit  tangent  to  it,  then  the  velocity  of 
propagation  of  edge  O'i  will  be  normal  to  VO*  as  well  as 

S>;  . 

•  nc  =o  •,  -  A-  =o 

The  propagation  will  be  determined  by  some  physical 
desirable  condition  of  energy  dissipation  in  the  growth  of 
the  crack,  which  may  be  deduced  by  considering  the  surface 
free  energy  or  cohesion  or  some  other  criterion.  This  is 
known  as  crack  propagation  criterion  and  will  be  discussed 
in  detail  in  the  next  section. 

Moreover,  in  accordance  with  the  causality  condition 
the  crack  propagation  velocity  should  not  be  greater  than 
the  P  wave  velocity, 

l£l  <irf 


Also  outside  the  crack  the  displacement  discontinuity  should 
van i sh , 


All-  -O 

^  2 

4.  1 

The 

surface  of  the  rupture 

should  be  free 

from 

traction; 

this  gives  the  following 

cond i t i on 

cTlT  '  n ;  -o 

C. 

4.2 
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Therefore  the  problem  is  specified  by  the  equation  of 
motion 


4  0 

^  LLl  —  (Ccj  jpty  5^-0 


4.3 


subject  to  the  boundary  conditions  4.1  and  4.2  together  with 
the  rupture  criterion  and  the  initial  conditions 


tie  =Uc  =0  t  <:o  4  4 

In  the  above  formulation,  to  solve  the  equation  of 
motion,  prior  Knowledge  of  and  the  normal  Vl[  to 

■2T(.X;  tj  are  required.  For  the  sake  of  simplicity  it  is 
assumed  that  the  crack  developes  along  some  preferred  plane 
so  is  known.  Also  assuming  that  the  growth  of  the 

crack  takes  place  in  some  specific  mode  and  with  knowledge 
of  the  propagation  velocity  \y  ,  the  dislocation  surface 
1)  may  be  defined. 


In  general  the  surface  may  be  irregularly 
shaped  but  the  restriction  to  the  crack  surface  imposed 
above  makes  the  solution  for  the  displacement  field  not  so 
appropriate  for  the  near  field  region,  where  the  details  of 
the  surface  may  have  profound  influence  on  the 
radiation  field.  However,  for  the  far  field  region  these 
minute  details  which  have  been  excluded  donot  have  any 
influence  on  the  displacement  field. 


s5 
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There  are  three  basic 
and  Eshelby,  1968),  in 

X3  =  0  and 


modes  of  the  crack  growth  (Bilbly 
each  the  crack  occupies  the  plane 

<$  -<>o^Xa<  do  (  fig.  4.1  ). 


T ens  i  1  e 


I  n-p  1 ane 


Ant i -p 1 ane 


Fig  4.1  Three  basic  modes  of  crack  propagation. 

I  Tension  mode 

The  state  of  tension  mode  of  crack 
opening  is  defined  by  on  the  crack  plane,  ie, 

X^-O  together  with  only  one  non  vanishing  component  of 
initial  stress  namely  ^  33  . 

I I  In-plane  mode 

The  faces  of  the  crack  slide  past 
one  another  as  indicated  by  the  arrow,  the  crack  opening  is 
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defined  by  &  Ul  ~  ■ f/(X,)  on  X3=0  and  the  only  non  zero 

initial  stress  is  (T,^  . 

Ill  Anti-plane  mode 

The  faces  of  the  crack  move  in  and 
out  of  the  figure, as  indicated  by  the  arrow  head  and  tail, 
here  AU_2  ~  ^  C^i)  on  X3=0  with  the  only  non  vanishing 
initial  stress  Cf£3  . 


Tension  crack  opening  mode  is  appropriate  for  regions 
under  tension  stress,  such  as  ridges,  but  at  a  depth  below 
few  hundred  meters  lithospheric  pressure  would  keep  the 
faces  together,  so  a  suitable  mode  of  rupture  may  be 
in-plane  or  anti-plane  shear  crack  or  in  general  a 
combination  of  both. 


The  problem 

assumes 

that  the  growth  of 

crack 

starts 

at 

0  .  Str ict ly 

speaking 

the  presence  of  a 

sma  1  1 

crack 

at 

the 

origin  at 

initial 

time  is  assumed, 

a  1  though 

the 

dimension  of  that 

sma  1  1 

opening  is  not  taken 

i  nto 

the 

account  in  mathematical  formulation. 

Further,  invoking  symmetry  and  anti  symmetry  of 
stresses  and  displacement  fields  about  the  plane  X3-0  , 
the  full  space  problem  may  be  reduced  to  a  b.v.p.  for  half 
space.  Thus  the  solution  has  to  be  found  in  only  one  region 
and  using  symmetry  arguments  the  solution  in  the  other 
region  can  also  be  found. 
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Most  of  the  work  done  so  far  in  seismic  source  dynamics 
consists  of  solving  the  problem  formulated  above  for 
different  growth  modes  with  different  crack  propagation 
cr i ter i a  . 


4  .  1  Crack  propaqat i on  cr i ter i a :  - 


To  describe  the  process  of  rupture  as  a  spontaneous 
phenomenon,  one  has  to  specify  the  change  in  dimension  of 
the  crack  in  addition  to  the  boundary  conditions  at  the 
surface  of  the  crack.  Griffith  (1920)  described  this  change 
in  dimension  by  considering  global  energy  balance  in  the 
material.  Barenblatt  (1962)  considered  a  yielding  mechanism 
at  the  crack  tip  and  developed  another  fracture  criterion. 

Griffith  energy  ba 1 ance  mode  1 : - 


Griffith  (1920)  formulated  crack  propagation  criteria 
in  terms  of  energy  balance  of  the  whole  material  surrounding 
the  crack.  There  are  always  weak  zones  inside  a  material  and 
under  the  applied  load,  if  the  rate  of  release  of  mechanical 
energy  (applied  load  work  done  W L  plus  strain  energy  JE 
stored  in  the  medium)  exceeds  a  certain  minimum  value, 
growth  of  these  zones  takes  place  resulting  in  the 
development  of  the  cracks  whereby  new  surfaces  having  some 


. 


surface  energy  \JS  are  created.  Therefore,  the  total  energy 
of  the  system  is, 
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U  -  (r  WL-+  U&)  +  U  4.5 

The  condition  for  a  crack  to  grow,  therefore  is  that 
the  decrease  in  mechanical  energy  should  be  just  greater 
than  the  increase  in  surface  energy.  Thus,  under  equilibrium 
condi t ion 


where 


d  i s  a  measure  of  crack  length. 


4.6 


Clearly,  the  crack  would  extend  or  close  up  reversibly 
for  small  displacements  from  the  equilibrium  length, 
according  to  whether  the  left-hand  side  of  4.6  is  negative 
or  pos i t i ve . 


The  theoretical  calculation  as  well  as  observed  facts 
suggest  that  there  is  stress  concentration  near  the  tip  of 
the  crack  and  the  growth  of  crack  occurs  when  it  exceeds  a 
certain  critical  value.  Irwin  (1958)  introduced  the  concept 
of  stress  intensity  factor  (the  stress  at  the  tip  of  an 
infinitely  sharp  crack  due  to  the  remote  stress),  K  ,  and 
extended  the  Griffith  concept  to  account  for  the  behaviour 
of  real  materials  by  considering  a  dissipative  component  in 
the  potential  energy  term  to  account  for  nonlinearity  at  the 
crack  tip. 


' 


One  defines  the  crack  extension  force  (the  mechanical 
strain  energy  release  rate  per  unit  width  of  the  crack 
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front)  as  (Lawn  and  Wilshaw,  1975,  p.50), 

dL(.-WL-rUE) 

^  ~  ~~  d  a 

which  is  related  to  K  through  elastic  constants  . 

One  also  defines  the  fracture  surface  energy,  P  ,as 
(Lawn  and  Wilshaw,  1975,  p.79), 

ir  =  P5 

d  Cl 

At  Griffith  equilibrium,  'the  so  called  Irwin  and  Orwan 
generalization  of  the  Griffith  criterion'  (Lawn  and  Wilshaw, 
1975,  p .  79 )  is 

Gc  =  “2.  r 


where  the  subscript  c  refers  to  the  equilibrium  value. 

So  fracture  will  extend  if 

If  there  is  no  dissipative  component  in  the  work  of 
creating  new  fracture  surface, 

r  =y 

where  Y  is  the  reversible  surface  energy  of  an  ideal 
brittle  solid. 

As  the  crack  starts  to  propagate,  material  near  the 
sides  of  the  crack  is  displaced.  The  rate  at  which  this 


i  m  .e  >n< 
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material  can  move  limits  the  velocity  of  the  crack 
propagation.  Following  the  Mott  (1948)  procedure,  a  Kinetic 
energy  term,  Li  ,  is  added  to  the  expression  for  the  total 
energy  and  invoking  the  equation  4.6  the  fracture  criterion 
i  s : 

G  -  >  If  4.7 

A  & 


Cohesive  zone  model 


Barenblatt  (1962)  developed  the  cohesive  zone  model  by 
considering  crack  propagation  on  the  atomic  level,  where  it 
may  be  visualized  as  breaking  of  the  bonds  between  the  atoms 
at  the  crack  tip.  The  force  holding  atoms  together  is 
plotted  against  inter-atomic  separation  (fig.  4.2  ). 


Pig.  4.2  Cohesive  force  between  atoms  as  a  function  of  their  separation  distance 
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At  equilibrium  separation  e  the  force  is  zero;  with 
the  increase  of  the  separation  it  grows  and  reaches  a 
maximum  for  the  separation  and  after  that  it  diminishes 
rapidly  with  further  increase  and  beyond  it  may  be 

regarded  as  zero.  Therefore,  it  may  be  assumed  that 
beyond  ^  a  weakening  of  bonds  starts  and  at  ^  it  is 
completely  broken.  Thus  at  the  tip  of  the  crack  the 
inter-atomic  separation  will  be  of  the  order  of  and 
cohesive  force  will  be  maximum  there,  and  trailing  behind 
will  be  a  zone  with  steadly  increasing  separation,  which 
will  run  into  the  completely  fractured  region  with  zero 
cohesive  force  beyond  the  point  where  separation  is  of  the 
order  of  (fig.  4.3  ). 


Fig.  4.3  Atomic  view  of  the  fracture  process. 


Clearly  the  gradient  of  the  cohesive  force,  ie,  stress  field 


will  have  a  hump  which  is  of  the  order  of  a  few  atomic 
distance  over  the  transition  zone.  The  amount  of  work  to 
produce  unit  area  of  crack  between  originally  unseparated 
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planes  of 


atoms  equal 


to 


<rc  =  n 


ft*)  s 


must  be  done  against  the  cohesive  force;  where  D  is 
number  of  bonds  per  unit  area.  With  this,  the  expression 

the  stress  intensity  factor  (Bilbly  and  Eshelby,  1968) is: 

o 


the 

for 


where  d  is  the  dimension  of  region  over  which  the  cohesive 
force  exists.  This  gives  a  finite  and  continuous  stress  at 
the  crack  tip. 


4 . 2  Dynami c  energy  ba 1 ance : - 


Using  the  overall  energy  rate  balance,  the  energy  flux 
into  the  crack  tip  of  an  extending  crack  can  be  calculated. 
The  analysis  here  follows  Atkinson  and  Eshelby  (1968)  and 
Freund  (1972).  For  simplicity,  let  us  consider  a  plane  crack 
(fig.  4.4)  in  a  2D  deformation  field,  extending  at  an 
arbitrary  rate  {IS  )  where  the  figure  represents  the  body  at 


■ 
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a  certain  fixed  instant  of  time. 


Suppose  that  the  crack  is  extending  through  the  motion 
of  either  or  both  crack  tips.  At  any  instant  of  time,  the 
rate  of  work  (  P  )  of  the  tractions  on  5  is  equal  to  the 
rate  of  increase  of  internal  energy  of  the  body  plus  rate  at 
which  energy  is  absorbed  by  the  moving  tip.  By  the  body,  is 
meant  the  region  V  in  the  limit  as  S,  and  shrink  on 
to  the  crack  tip.  Denoting  the  total  crack  tip  energy  flux 
by  ot  +  iv  the  energy  rate  balance  is 

?  ^  UE  4-^ 


Taking  into  the  account  that  5,  and  move  with 
the  crack  tip  and  expressing  each  term  of  the  above  equation 
in  terms  of  elastic  fields,  we  get: 


where  Y\-  is  the  normal  to  %.  directed  away  from  the 
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crack  tip. 


Clearly  the  energy  flux  depends  only  on  the  near  tip 
elastic  field  and  the  calculated  value  of  is 
independent  of  choice  of  contour  (Freund  (1972)  ,eq. 


Xv 


14).  Let  (Ut),  t )  be  the  position  of  the  crack  tip  at  any 
instant  t,  and  choose  5*  as  a  rectangular  segment 

I X,  —  £&)  1  =  £>,  and  |  Xal—  &2_  *  where  6t  and  are 

abritrarly  small  positive  numbers.  If  the  limit  as 
is  interpreted  as  the  ordered  limit  °  followed  by 

£.->0  ,  then  the  only  contribution  to  the  integral  in  the 

equation  4.8  comes  from  the  segments  parallel  to  the 
X-axis,  on  which  L %\~0  and  the  equation  4.8  reduces  to 


& 


Lt  (  t(x/,o)|j?  <f)  -  u(x',  ,CT)]  elk/ 


4 


UrW.  x/  -  X,  —  {ft 


4.9 


The  energy  flux  (g)  is  related  to  the  crack  extention 


force  (G)  by, 


3  =  2  0-G, 


where  the  factor  2  accounts  for  both  faces  of  the  crack. 
Using  the  equation  4.9  we  find  (  Achenbach,  1974  ) 
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c  -  JS  Altf)  BO*) 

A  4.10 

where  A(v)  B ( v )  are  the  functions  only  of  the  crack 
propagation  velocity  that  decrease  from  a  value  near  1  at 
ir=0  to  zero  at  the  terminal  velocity.  For  anti-plane  cracks 

v  -'4 

Alfr)  =  0  -  ;  Bdr)  =  C  1  + 

and  for  plane  cracks 

to)- s(-Xr)jpr^  >  _2(47i£-'J 


I 

Here  is  the  Rayleigh  wave  velocity,  and  S ( p )  is 

the  de  Hoop  (1958)  function  defined  as 


500  -  £kj=|-i  |  ^ 


--r.) 


and  the  stress  intensity  K  may  be  interpreted  as  the 
stress  concentration  that  would  remain  if  the  crack  suddenly 
stopped;  it  is  a  function  only  of  the  stress  and  extension 
history  of  the  crack  but  it  is  independent  of  IT  (  Fossum  and 
Freund,  1975;  Freund,  1976).  As  G  is  related  to  material 
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constant  (Irwin  criterion)  so  the  equation  4.10  is  also  the 
equation  of  rupture  velocity.  Kostrov  (1970)  has  developed 
following  equation  governing  rupture  velocity 


where 


and  subscript  to  the  stress  intensity  factor  term  refers  to 
the  corresponding  rupture  mode.  Here,  a  general  situation 
involving  all  three  modes  of  rupture  is  considered.  Freund 
(1979)  has  given  an  extensive  review  on  the  mechanics  of 
dynamic  crack  propagation. 

4 . 3  Shear  Cracks : - 

Perhaps  the  simplest  models  for  earthquakes  are 
transient  shear  cracks.  At  time  t=0  suddenly  a  crack  appears 
at  the  origin  and  starts  growing  obeying  the  equation  of 
motion 


(r*v  ^ 


i  til.  _  ^ 

-  - —  7  -L.  '  n 


4.  12 


for  anti -plane  crack  and 


4.13 


. 
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for  plane  crack,  with  the  boundary  conditions 


T1J=-  KM)  =  «r£5-<r£  ;X3=0 


4.14 


U.S.  =  0 


X,<(- ,  v,  >1  +  -j  *3  =o 


for  antiplane  and 


X,3  r  -  KVt)  =  -  ^,3 


L(t;  cx,  <Kfc)  j  ^3 


4.15 


u.,  =  0 


X,  <1,  ,  X,  >£+  x5.o 


for  plane  crack,  to-gether  with  the  initial  conditions 


u.z  ^  az  -  o 


t  <  0 


4.  16 


for  antiplane 


u, -a,  -  a3^  u3-  0 


4.17 


for  plane  case. 

Due  to  the  symmetry  w.r.t.  the  plane  0  it  is 

sufficient  to  find  the  solution  only  in  one  half-space. 
Hereafter  the  displacement  on  the  crack  plane  will  always 
refer  to  that  determined  for  the  half-space.  The 
displacement  discontinuity  (slip)  will  be  twice  the  surface 
displacement  determined  for  the  half-space. 

Using  Green's  function  technique,  the  integral  equation 
for  displacement  is  given  by  (  Burridge,  1969  ) 


4.18 


. 


■ 


'  ■ 
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For  anti -plane  crack  the  Green's  faction  for  the 
half-space  <0  is  given  by  (  Achenbach,  1974  ) 


H  [Ct'-t)  -  Cx/-*tj / us 

^2.1  ,  cl.  r  .  ,  -,2- 1  V, 


TTf'CCt'-t) 


A 


4.  19 


With  that,  the  solution  is 


U  rL\  , 


/  A  ,  / 


l 

tr^) 


/Y 


Tx3CXi,  ^X|  ^ 


4.20 


where  region  7  of  integration  is  the  backwards 
character i st i c  triangle 

<  &  (t-t) 


Fig.  4.6 


. 
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Region  f>0  may  contain  an  area  outside  crack  where  ^13 
is  not  known.  So,  to  solve  above  equation,  TX3  has  to  be 
found  in  that  region.  Kostrov  (1966)  has  given  it  in  terms 
of  the  stress  drop  inside  the  crack  assuming  that  a 
contribution  comes  only  from  one  of  the  crack  tips. 


4.21 


where  L(t')  is  the  location  of  the  crack  tip  and  “fc,.  is  the 
time  when  the  crack  tip  locus  intersects  the  integration 
path,  namely 

Uk^-x,'  -  ur4ta  -  i(W 


Considering  the  energetics  at  the  crack  tip,  Kostrov 
gave  the  following  equation  governing  the  location  of  the 
crack  tip 

l 


JT- 


X, 


=  (2  7T/U  6] 


V 


2- 


VI  -  w/tr* 


4.22 


The  above  equation  holds  only  when 

€ 

,  &*•  ' 


k 


2- 


4.23 


Using  transformation 


■ 
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The  equation  of  the  motion  for  the  crack  tip  can  be  written 
as 


CUT/M&,) 


/7_ 


-u-x, 

t-  ust 


4.25 


Once  the  locus  of  the  crack  tip  is  determined  T13 
can  be  calculated  and  the  displacement  field  on  O  is, 

using  a  transformation  similar  to  4.24,  given  by 


/ 


4.26 


Thus  the  displacement  discontinuity  or  slip  will  be(-2Uj  and 
the  slip  velocity  (  Achenbach,  1974  ) 


to 


AlUx„t)  = 


L r 


<1 


4 


-X,)2- 


Kx,,  tr- (£ t) — x, 

flit)  - 


Ax, 


o 


4.27 


Madariaga  (1977),  using  the 


above  equations,  did 

<9-51^  (fig. 


calculations  for  the  elastic  field  with 


' 


56 


4.7). 


vst 


Static  field  just  before 
initiation 


Stress  change  due  to  rupture 


Total  stress  on  the  fault  plane 


Slip  velocity  distribution 


Fig.  4.7  Elastic  field  on  the  plane  of  suddenly  starting 
anti-plane  crack  (Madariaga,  1977). 


For  in-plane  we  have  the  solution 


5. 


Gnn  (X  At 


4.28 
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and  the  Green's  functions  given  by  Lamb  (1904)  are 
fol lows . 


as 


(/-  tr;  V- 


TT/UA^X, 


X  i 


R  (X) 


-L  <X  < 
LXC 


& 


_1 

U, 


<  cX 


4.29 


G„(\,o,-t)  =  —  SCt  - 
3  ^ 


___  e^-ir 

TTpr^X, 


-+ 


&%) 


where  <  -  £  ,  and  R  is  the  Rayleigh  function 

RfeU  -  (<*-  Off* 

fUO  -  (5^-l h)  +Aa^ (■<’’- 


and  KL is  defined  by 


^  , the  area  of  the  integration,  is  backwards 
character i st i c  triangle 

x/-  X)J  <,  Up 

which  has  region  lying  outside  the  crack  for  which  Tj3  is 


58 


not  Known.  Similar  to  the  anti-plane  crack  Kostrov  (1975) 
obtained  the  analytic  solution  using  integral  transforms  and 
Wiener-Hopf  technique.  The  results  are  much  more  complicated 
than  the  anti-plane  case  and  valid  only  for  . 
fVladar  i  aga' s  (  1  977  )  calculations  based  on  those  results  are 
reproduced  here  (  fig.  4.8) 


Fig.  4.8  Elastic  field  on  the  plane  of  suddenly  starting 


in-plane  crack  (Madariaga,  1977). 
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Slip  dependent  boundary  condi t ion : - 


A  gradual  creep  deformation  prior  to  the  sudden  motion 
on  the  fault  is  observed  in  the  field  as  well  as  in  the 
laboratory  experiments.  Ida  (1973)  suggested  a  constitutive 
relationship  between  stress  and  slip  rate  to  be  incorporated 
in  the  solution  of  the  crack  propagation  to  take  into  the 
account  creep  motion.  He  gave  an  analysis  only  for  the 
anti-plane  case. 


Fig.  4.9  Constitutive  relation  between  stress  and  slip  rate 
(Ida,  1973).  is  critical  slip  rate 


<ac 

i  4.32 

As  tectonic  stress  increases  over  dynamic  frictional 
stress  (  <r^  )  creep  motion  starts  along  the  fault  plane  and 
the  rate  of  slip  increases  with  the  stress  increase  till  it 
reaches  a  critical  (static  frictional  stress)  value,  beyond 
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that  stress  drops  to  the  level  (T^  creating  an  earthquake. 

Let  us  consider  that  a  semi - i nf i ni te  crack  suddenly 
appears  for  at  t=0  .  Thus,  there  will  be  two  time 

intervals  and  •  F°r  t<--^  there  will  not  be  any 

tfs,  Uf  \r4 

disturbance  from  the  crack  tip  so  during  that  period  the 
slip  will  be  uniform  and  the  stress  drop  will  be  determined 
by  the  constitutive  relation  4.32. 


The  displacement  field  is  obtained  by  integrating 
(4.26)  over  appropriate  limits  and  yields 


a 

/X 


Using  relation  4.32 


r '  = 


C  +  Xx: 


=  <r. 


o 


Thus  the  stress  drop 
material  constant  is 


(=  -  r ai1=  - 2. r  1,  - 

in  terms  of  initial  stress  and 


i + r  Ifi 


4.33 


For  t  we  have  disturbances  coming  from  the  crack  tip, 
so  both  )?  as  well  as  are  unknown  and  the  solution  can 
be  obtained  by  solving  both  4.26  and  4.32  simultaneously.  To 
solve  it  Ida  (1973)  developed  a  numerical  scheme  based  on 
discretization  of  integral  equation. 
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j 

For  various  choice  of  parameters  (J£  ,  f  and  (f^  Ida 
made  numerical  calculation  and  found  that  there  are  two 
distinct  classes  of  rupture;  in  one  of  them  rupture  once 
started,  accelerates  smoothly  to  \%  (fig.  4.10)  and  in  the 
another  class  rupture  after  propagating  comes  to  rest,  then 
restarts,  stops,  and  repeats  the  process  (fig.  4.11). 


Fig.  4.10  Smooth  rupture  propagation  r-  ■  _  .  ..  ,  , 

K  H  y  Fig  4.11  Irregular  rupture  propagation 


(Ida.  1973  )  . 


(Ida,  1973) 


Distribution  of  (TT  in  the  unit  of  F5  Time  is  measured  in 

3  Jo 

the  unit  of  tD  and  distance  in  the  unit  of  l^tD  Where 

is  the  time  required  for  crack  initiation  CTT  Gi  / 1-  \J±  )  . 


Andrews  (  1976)  ,  for  in-plane  cracks  simi lar  to  Ida  s, 
developed  a  slip  weakening  model  assuming  that  the  shear 
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stress  across  the  crack  plane  is  a  function  of  the  slip  and 
when  the  stress  reaches  the  static  friction  limit, 
instability  begins  and  weakening  occurs.  For  slip  greater 
than  ci  ,  the  stress  drops  to  the  dynamic  friction  and  the 
work  done  up  to  that  point  is  identified  as  specific 
fracture  energy, 


With  this  model  he  examined  the  rupture  propagation  by 
a  finite  difference  scheme.  The  results  are  dicussed  in 
terms  of  two  nondimensiona  1  numbers,  Lc/[_  and 

5  -  - Z.-—  ,  where  o-0  is  the  initial  stress  and  I  is 

<r°  _<r^ 

the  critical  half  length  of  the  in-plane  Griffith  crack 


which  is  given  by 


v°V 
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The  calculation  shows  that  when  the  parameter  S  is 
greater  than  1.63,  then  is  always  less  than  , 
and  ^  approaches  as  the  crack  length  increases. 


Fig.  4.13  Rupture  propagation  with  slip  weakening  model 
( Andrews ,  1976), 

But  for  S  =  1 . 63  the  crack  is  accelerated  towards  and  may 
approach  Op  as  the  crack  length  increases.  The  range 
appears  to  be  forbidden.  Burridge  (1979)  has  also 
confirmed  this  forbidden  range. 

Stopping  mechani sm  of  the  cracks : - 


Husseini  (1975)  suggested  that  along  trajectory  of  a 
running  crack  either  an  increase  in  fracture  energy  or 
finiteness  of  available  strain  energy  will  arrest  the 
propagation  of  the  crack.  An  example  of  the  former  situation 
may  be  a  rupture  initiating  on  pre  existing  fault  locked  by 
frictional  stress  and  attempts  to  propagate  into  the  region 
of  fresh  unfractured  rock  possesing  greater  fracture  energy. 
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And  a  region  in  the  path  incapable  to  supply  any  stress  drop 
to  the  crack  when  the  crack-tip  passes  them,  might  be  an 
example  to  the  latter  case.  For  semi  infinite  anti-plane 
crack  Husseini  examined  these  mechanisms  and  found  that  if 
arrest  is  through  the  presence  of  fracture  energy  barriers 
then  the  inequality  to  be  satisfied  is 


4.35 


where  AT  is  the  increase  in  fracture  energy,  M'is  the 
average  stress  drop  and  R  is  the  character i st i c  radius  of 
the  fault. 

If  on  the  other  hand  the  earthquake  is  arrested 
according  to  the  finiteness  of  strain  energy  or  seismic  gap 
model  then 


4.36 


Thus  if  arrest  is  predominantly  caused  by  barriers  then 
the  the  inequality  4.35  will  constitute  a  weak  relationship 
between  the  fracture  energy,  stress  drop  and  radius.  If,  on 
the  other  hand,  arrest  is  caused  by  seismic  gaps,  equation 
4.36  constitutes  a  strong  relationship  between  fracture 
energy,  stress  drop  and  radius.  Interestingly  both  4.35  and 
4.36  predict  that  the  fracture  energy  will  be  proportional 
to  the  radius  and  the  square  of  the  stress  drop. 
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The  equation  4.36  leads  to  the  following  relationship, 


The  combination  of  this  relationship  with  the  Keilis-BoroK 
(1959)  relationship  between  stress  drop  and  seismic  moment 
(  M  )  yields , 


This  relationship  is  statistically  acceptable  for  the 
Southern  California  faults  and  the  Tonga-Kermadec  Arc 
earthquakes  suggesting  that  the  seismic  gap  model  of  arrest 
may  be  valid  in  those  regions. 

Andrews  (1975)  applied  numerical  methods  to  the  two 
dimensional  shear  cracks  problem  and  found  that  with  an 
appropriate  nonuniform  initial  shear  stress  the  crack  would 
stop  by  i tsel f . 

Radiation  field  of  the  crack  mode  1 : - 

Richards  (1973,76)  solved  the  problem  of  a  3D 
prestressed  medium  in  which  the  rupture  front  starts  at  a 
point  and  propagates  out  as  an  ellipse  over  the  fault 
surface  at  constant  velocity.  From  the  considerat ion  of  the 
stress  near  the  fault  surface,  he  concluded  that  the  rupture 
velocity  will  be  between  (J^  and  t ^  .  The  synthetic 


seismograms  calculated  by  him  are  reproduced  below  (-fig. 
4.14). 
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Fig  4.14  Synthetic  seismograms  due  to  elliptical  rupture 
( R i chards ,  1976) 

Madariaga  (1976)  also  applied  numerical  methods  to  a 
plane  circular  crack  in  3D  which  expands  at  a  prescribed 
constant  velocity  and  then  suddenly  stops.  The  solution  for 
slip  on  the  fault  and  far  field  body  waves  were  obtained 
(fig.  4.15).  The  far  field  displacements  were  strongly 
affected  by  the  energy  radiated  when  the  fault  stopped. 
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I0~2  I0~'  I  10 

FREQUENCY 

Fig  4.15  Far  field  spectra  for  P  and  5  waves  due  to 
circular  crack  (Madariaga,  1976). 

0  is  the  take-off  angle  (angle  between  X-¬ 
axis  and  direction  of  observation). 


On  further  investigation  Madariaga  (1977)  found  that 
due  to  the  sudden  starting  or  stopping,  high  frequencies  are 
radiated  from  the  fault  and  are  determined  by  the  change  in 
the  slip  velocity  field.  Since  the  highest  slip  velocities 
occur  right  behind  the  rupture  front,  it  is  expected  that 
these  radiations  come  from  the  region  in  the  neighbourhood 
of  the  crack-tip.  The  radiation  pattern  due  to  the  sudden 
starting  or  stoping  of  a  shear  crack  is  shown  in  the  figure 
4.16.  Ak i ( 1980)  has  presented  an  extensive  review  on  dynamic 
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source  theory. 


stopping  in-plane  crack  (  P  and  SV  waves)  and 
anti-plane  crack  ( SH  wave)  (Madariaga,  1977). 

It  may  be  pointed  out  that  the  validity  of  the  result 
presented  here  is  much  restricted  for  a  real  earthquake  due 
to  the  simplifying  assumptions  considered  .  The  assumptions 
are 

1.  an  unbounded  solid 

2.  homogeneous,  isotropic,  linear  elastic  material 

3.  a  plane  two  dimensional  geometry 

4.  extension  of  crack  in  its  own  plane 

5.  constant  value  of  F 


6. 


brittle  fracture 
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Obviously  it  will  be  desirable  to  remove  some  of  the 
restricting  assumptions.  In  particular,  since  it  has  been 
observed  in  the  field  that  faults  may  bifurcate,  an 
extension  of  the  analytical  work  to  include  branching  and 
asymmetric  crack  propagation  deserve  primary  attention. 
Extension  of  analytical  work  to  include  interactions  between 
the  tip  of  the  crack  is  also  of  interest.  Also  the  plastic 
flow  in  conjuction  with  dynamic  fracture  and  possible 
functional  dependence  of  P  on  rupture  velocity  should  also 
be  the  subject  of  study. 


5.  SIMULATION  OF  FAULT  MOTION  AND  EARTHQUAKE  OCCURRENCE 


There  have  been  attempts  to  simulate  the  fault  motion 
and  the  occurrence  of  earthquakes  by  experimental  as  well  as 
by  numerical  means.  Generally  it  is  accepted  that  faulting 
associated  with  earthquakes  is  related  to  the  relative 
sliding  of  two  rock  surfaces  with  frictional  resistance 
opposing  the  motion.  In  the  simulations,  the  fault  regions 
are  repesented  by  discrete  blocks  which  are  driven  to  slide 
on  the  friction  surface  by  various  combinations  of  elastic 
and  viscous  force  elements. 

In  the  experiments  on  rock  friction  it  has  been  found 
that  there  are  two  modes  of  displacement  along  a  frictional 
surface:  stable,  in  which  there  is  continuous  displacement 
without  abrupt  change  in  the  applied  force,  and  stick  slip 
which  is  a  jerky  motion.  Under  the  varying  condition  of 
normal  stress,  displacement  rate,  pressure,  temperature, 
pore  fluid  etc.  both  occur.  Also  transition  from  one  mode  to 
another  takes  place  but  this  is  not  usually  abrupt.  At  low 
pressure,  say  10  bars,  the  sliding  mode  is  seen  to  be  stable 
and  changes  to  stick  slip  at  higher  pressure.  With  increase 
of  temperature,  due  to  ductile  flow  the  transition  from 
stable  to  stick-slip  is  inhibited.  It  has  been  observed  in 
sandstones  that  with  the  increase  of  velocity  a  transition 
from  stick-slip  to  stable  takes  place.  Also  higher  pore 
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fluid  pressure  promotes  stable  sliding  (Logan,  1975). 

The  motions  on  the  fault  surface  are  found  to  take 
place  suddenly.  During  the  sudden  slip,  shear  stresses  are 
relieved  and  the  fault  surface  may  thus  remain  locked 
together  until  some  time  later  slip  suddenly  takes  place 
again.  Such  sudden  intermittent  motion  on  a  pre-existing 
fault  in  the  earth  is  similar  to  jerky  sliding  between  the 
rock  surface  in  laboratory  experiments  (Brace  &  Byerlee, 
1966).  Therefore  stick  slip  is  consider  to  be  a  possible 
mechanism  of  slip  along  the  fault.  The  assumption  made  to 
explain  stick-slip  is  the  fact  that  dynamic  friction  is  less 
than  static  friction.  The  physical  basis  put  forward  for 
this  are:  thermal  softening,  creep  instability  and  brittle 
failure.  Byerlee  (1970)  examined  these  alternatives  in 
detail  and  suggested  that  the  most  probable  cause  may  be 
brittle  fail  ure . 

It  has  been  observed  that  there  is  an  increase  of 
temperature  at  the  sliding  surfaces.  Since  a  rise  of 
temperature  lowers  the  strength  of  material,  this  leads  to 
the  idea  that  friction  during  s 1 idi ng , i . e . ,  dynamic  friction 
may  be  less  than  static  friction.  If  this  is  the  only  cause, 
then  before  stick  slip  starts  there  should  be  always  a 
period  of  stable  slip  during  which  temperature  rises  to 
decrease  the  strength  of  material.  This  is  not  true. 

When  surfaces  are  held  to-gether  they  make  contact  only 
at  a  few  points.  The  creep  instability  model  assumes  that 


. 
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these  points  of  contact  or  junctions  deform  by  creep 
mechanism  so  that  the  size  of  the  junctions  increase  with 
time.  During  the  sliding,  the  time  of  contact  is  so  small 
that  junction  growth  does  not  occur.  If  the  friction  force 
is  the  force  required  to  shear  the  junctions  then  the  force 
required  to  start  sliding  after  the  block  has  been  at  rest 
may  be  greater  than  the  force  required  to  maintain  sliding. 
Hence  static  friction  could  be  greater  than  dynamic 
friction.  Dieterich  (1972a)  has  observed  this  effect  and 
cone  1 udes 

I U  (£ )  jU  (p )  \^  \  -f  c(  Urt^  ^  tj  At  >  I  £e  Ccrv'i  5  .  1 


where  oC  i s  a  constant  and  is  the  duration  of  the  contact 
following  sliding.  Recently  on  the  basis  of  the  laboratory 
data  he  (Dieterich,  1979)  has  developed  the  following 
quantitative  relationship  between  cofficent  of  friction  /a 

9 

and  time  of  stationary  contact  and  slip  velocity  li(t)  . 


-+■ 


a] 


5.2 


where  C,3Cx,Ca  ^  ^  ^  are  constants. 

The  brittle  instability  model  for  stick-slip  considers 
that  asperities  along  the  fault  surface  lock  together  to 
form  the  frictional  resistance.  When  the  shear  force  becomes 
strong  enough  to  cause  brittle  failure  of  these  asperities, 
motion  can  occur.  As  the  sliding  proceeds,  continual  locking 
and  breaking  of  asperities  causes  fluctuations  in  the 
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frictional  resistance.  Thus  the  frictional  resistance  varies 
in  some  irregular  manner  with  displacement.  The  value  of 
frictional  resistance  at  the  initiation  of  an  unstable 
sliding  event  is  analogous  to  the  static  friction,  and  the 
average  value  of  fluctuating  friction  during  the  sliding 
corresponds  to  dynamic  friction  (Nur  and  Shultz  (1973),  Nur 
(1978)).  The  presence  of  gouge  along  a  sliding  surface  is 
strong  evidence  in  the  favour  of  this  brittle  instability 
mechanism  of  stick-slip. 

The  modeling  of  fault  motion  also  requires  some 
relationship  between  the  stress  acting  on  the  rocks  and 
the  consequential  strain  .  In  general  this  relationship  can 
be  quite  complicated,  but  even  by  restricting  the  analysis 
to  a  linear  relation  between  them,  considerable 
understanding  of  the  process  can  be  had.  Thus  for  a 
connecting  element  of  the  model,  the  stress-strai n 
relationship  can  be  written  as 

i r-Ofr 

# 

where  0  is  a  linear  operator  on  £  . 

For  an  elastic  spring,  0  &  =•  dbt  =  (<£  ,whi  le  for  a 
vi scous  element ,  Oe  ' ^ 6  (Jaeger  and  Cook,  1976). 

More  complicated  connecting  elements  can  be  constructed 
by  the  combination  of  these  simple  elements,  and  applying 
linear  operator  techniques,  stress-strain  relationship  can  be 
derived  (fig.  5.1a).  The  time  dependent  behaviour  of  stress 
and  strain  for  those  substances  is  shown  in  the  fig.  5.1b. 
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Fig.  5.1a  Multiple-element 
rheologic  model  (Cohen.  1979) 

(a)  Maxwell  substance 

(b)  Kelvin  substance 

(c)  Three  element  generalized 
Ke Ivin 

(d)  Alternative  three  element 
substance 

(e)  Four  element  Burgers  substance. 


Maxwell  substance 


Fig.  5.1b 

Time  dependent  behaviour 
of  stress  and  strain  for 
selected  linear  substance 


(Cohen,  1979). 


Generalized  Kelvin  substance 
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5 . 1  Simulation  mode  1 s  for  earthquakes : - 


Different  authors  have  tried  to  stimulate  stick-slip  by 
the  help  of  mechanical  models.  These  models  can  be 
categorized  into  two  classes:  the  first  consists  of  a  mass, 
sliding  over  a  frictional  surface,  driven  by  a  constant 
force  applied  on  the  block  through  a  spring  (Rider  model), 
another  consists  of  a  linear  array  of  masses  conected  via 
spring  elements  and  driven  by  springs  conected  to  a  moving 
plate  (Burridge  and  Knopoff  model).  Assuming  different 
friction  laws  and  rheological  models  for  connecting 
elements,  the  equations  of  motion  of  the  blocks  are 
constructed  and  then  rupture  velocity,  time  of  slide, 
frequency  of  slide  etc.  are  calculated.  These  results  are 
then  matched  with  natural  events.  With  these  some  authors 
have  also  constructed  laboratory  models.  A  systematic  detail 
of  different  models  is  presented  below. 

R i der  mode  1 : - 


Rider  model  consists  of  a  mass  spring  system  as  shown 
in  (fig.  5.2).  A  block  of  mass  YC\  sitting  on  a  rigid 
frictional  surface  is  subjected  to  stress  due  to 
accumulation  of  elastic  strain  within  a  spring  which  is 
being  stretched  at  constant  velocity  V  .  In  absence  of  any 
inertial  motion  in  the  block,  the  friction  force  prevents 


■ 
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any  sliding  until  the  stress  rises  to  the  maximum 

r5 

strength  j-  .  Once  sliding  has  begun,  the  friction  resistance 

ol 

drop  to  the  dynamic  value ,  which  tend  to  accelerate  the 
block.  This  results  in  a  decrease  of  force  in  the  spring  and 


1C  V 
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Fig.  5.2  Rider  model. 

the  block  eventually  stops.  Now  the  block  remains  stationary 
until  the  forces  in  the  spring  once  again  build-up  to  a 
value  to  overcome  friction  forces. 


Assuming  the  initial  position  of  the  block 
X-0  ,  motion  is  initiated  at  a  time  tD  when 

Ki/ 10  —  -f 


is  at 


5.3 


and  then  block  starts  moving;  the  equation  of  motion  of  the 
block  can  be  written  as 


wi  [i  =  v<  (yt  -  u/)  -  \ 


5.4 


the  solution  of  above  equation  is 


nit)  =  ji-  i/(^)  5^ (t'Q 
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For  practical  purposes,  once  the  sliding  starts,  the  block 
begins  to  move  much  faster  than  1/  ,  so  we  can  neglect  the 
effect  of  V  in  the  equation  5.4  .  This  is  equivalent  to 
setting  t'tc  in  the  equation  5.4  .  Thus  we  have, 

vna  -+-K ul  =  Kvt0  -  ~  -f5-  ^  M 


and  the  solution  is 


u.  ~  ~  0  -  ^  ( 17  ' 


so  we  have  velocity  of  the  block 


5.5 


a  - 


.=  ±1  s*(K)V-t. 

\fv<  iv\ 

and  the  acceleration 

=  X  UofeiXt't 


5.6 


» • 

IX 


wv. 


5.7 


Duration  of  slip  is  obtained  by  setting  in  5.6 

-  A  b  -  TT  /w/K 

Total  displacement 

__ 


AX 

Force  drop 


K 


A  F  ~  AlA'K-  —  "2-  ^ 

Peak  ve loci ty 


ti 

Peak  acceleration 

and  Potential  Energy  drop 


d. 


AE  =  I  AO- 
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Although  this  model  is  an  oversimplification  of 
complicated  tectonic  and  geologic  processes,  it  is  still 
able  to  give  a  qualitative  picture  of  the  fault  motion. 

Nur  (1978)  introduced  non  uniform  friction  in  the 
equation  of  motion  of  the  rider  model  to  incorporate  the 
brittle  instability  mechanism.  Since  in  the  rider  model 
there  is  no  distiction  between  displacement  or  position, 
therefore  let  J(x)  be  a  position  dependent  function 
representing  non  uniform  friction.  Then  the  equation  of 
motion  of  the  rider  can  be  written  as 


m  lb  -hKUl  -f(x) 


with  initial  conditions 


U.(o)  -  x„  ,  iU°)  =  V 


and  we  get  expression  for  time  of  the  rider  at  any  position 


On  inverting  above  equation  one  gets  the  displacement  time 
function  LUt')  . 

In  the  derivation  of  the  above  equation  the  nature  of 
non-uniform  friction  is  not  specified.  Thus  depending  upon 


the  nature  of 


a  solution  will  be  obtained  either  by 


analytical  or  numerical  means. 


Whitehead  &  Gans  (1974)  have  examined  the  behaviour  of 
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the  rider  model  with  the  assumption  that  the  rider  and 
sliding  surface  are  separated  by  a  thin  temperature 
dependent  viscous  layer  so  that  the  net  friction  forces 
decrease  with  increasing  velocity  and  the  solution  obtained 
shows  violent  motion  over  a  small  period  of  time  followed  by 
a  quiet  period. 

The  most  detailed  quantitative  calculation  on  rider 
model  has  been  done  with  variable  friction  and  a 
viscoelastic  element  by  Cohen  (1978).  The  introduction  of 
this  element  permits  modeling  of  the  transient  anelastic 
deformation  in  response  to  stress  loading  and  relaxation, 
and  provides  a  mechanism  for  partial  stress  recovery 
following  an  earthquake.  As  a  consequence,  it  is  possible  to 
simulate  the  episode  of  stable  sliding,  tertiary  creep 
preceeding  earthquake  and  long  term  aseismic  creep  also. 

Bur r i dqe  and  Knopof f  Mode  1 : - 


There  has  been  another  independent  effort  to  formulate 
dynamical  equation  for  particle  displacement  on  a  fault  by 
Burridge  and  Knopoff  (1967).  This  consists  of  a  number  of 
masses  connected  together  by  springs  and  resting  on  a 
frictional  surface.  The  masses  are  driven  by  springs 
connected  to  a  moving  plate  (fig.  5.2).  The  various  springs 
corresponds  to  elastic  and  viscous  properties  of  the 
material.  The  driving  plate  is  moving  with  velocity  V 
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which  is  unaffected  by  the  frictional  resistance  along  the 

fault.  Let  Fj  be  the  elastic  force  acting  on  the  i th 

f  5 

block  which  is  associated  with  static  friction  y. 

c 


Fig.  5.3  One  dimensional  Burridge  and  Knopoff  model 


Initially  all  blocks  are  at  rest  and  F:  <  ^ •  for  each 
block.  As  time  advances,  F '  grows  due  to  the  motion  of  the 
driving  block  which  increases  the  driving  spring  tension. 
Once  F-  overcomes  the  frictional  resistance  the  motion  of 
block  i  starts.  This  is  governed  by  the  following  equation 


ivi;  ,  F-  i; 

de 


J. 


— Fc  (u;+,  -  uu) -U;)  -+kc'  -  j-/ 


ft  ft 


Here 

Wl  =  mass  of  i'  th  block 

Kl*  =  spring  constant  for  the  spring  connecting  blocks  i 
and  i+1 

=  spring  constant  for  the  spring  connecting  i'th  blocks 
and  the  driving  plate 

•[j  =  dynamic  friction  resistance  of  the  i'th  block 

ii'  =  displacement  of  the  i'th  block 
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As  the  block  slides,  it  compresses  or  expands  the 
connecting  springs  to  the  adjacent  blocks  and  possibly 
stimulates  the  motion  of  these  blocks.  In  this  manner  slip 
will  propagate  along  the  fault.  Thus  the  motion  is  governed 
by  a  set  of  coupled  dynamic  equations  which  can  be  solved  by 
numerical  techniques.  Taking  different  numbers  of  blocks, 
different  values  of  masses,  different  rheological  elements 
and  friction  laws,  many  authors  have  examined  this  model 
(Burridge  and  Knopoff  (1967),  King  and  Knopoff  (1968), 
K i ng ( 1 975  )  , Ot suka  (1972),  Dieterich  (1972b, 73),  Cohen 
(1977),  Rundle  and  Jackson  (1977)). 

The  one  dimensional  model  for  laboratory  studies  of 
Burridge  and  Knopoff  (1967)  consists  of  eight  mass  connected 
by  a  spring  and  driven  by  stretching  of  the  spring  connected 
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Fig  5.4  Schematic  diagram  of  the  laboratory  model. 

to(fig.  5.4)the  first  mass. With  this,  it  is  found  that  small 
events  occur  largely  at  random  while  large  events  involving 
major  changes  in  the  elastic  potential  energy  are  almost 
periodic.  Between  major  events,  the  potential  energy  of  the 
system  increases  nearly  linearly  with  times  when  the 
intermass  springs  have  equal  spring  constants.  However  in 
case  of  unequal  spring  constants  this  is  not  so. 
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Fig.  5.5  Potential  energy  as  a  function  of  time  for  Burridge 
and  Knopoff  model  with  all  spring  equal  (Burridge 
and  Knopoff ,  1 967  )  . 


Furthermore,  the  major  events  occur  at  appoximately  the  same 
level  of  the  potential  energy  and  involve  comparable  energy 
drops  (fig.  5.5). 


The  plot  of  the  frequency  of  the  events  as  a  function 


Fig  5.6 

Frequency-energy  relationship 
with  all  spring  equal 
(Burridge  and  Knopoff,  1967). 


of  energy  released  (fig.  5.6)  except  for  the  lowest  energies 
fits  the  equation,  ,  _  ,  .  r 

L^10N=  cx-FU^fc 
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where  N  is  the  frequency  of  an  event  occurring  with  energy 
release  greater  than  or  equal  to  E  .  Assuming 

.where  MSl>1  -j  s  defined  as  magnitude  of  the 
simulated  events,  we  have 

L  M  -  a  -  b 

o  10 

This  is  analogous  to  the  Gutenberg-Richter  frequency 

magnitude  relationship. 

Typically  in  naturally  occurring  shocks  b  =  0.4, 
although  there  is  considerable  scatter  in  the  value  of  b 
But  in  the  above  model  b  =  1  Burridge  and  Knopoff 

attributed  this  discrepancy  to  the  one  dimensional  nature  of 
the  model.  Furthermore  the  natural  shocks  follows  a  Poisson 
distribution  but  in  this  case  some  nonPoisson  component  is 
present,  presumably  due  to  the  interaction  between  adjacent 
b 1 ocks . 

i  i  '  DRIVING  BLOCK  i 

,  ' 


A-FAULT  \  VISCOUS  '  P-FAULT 

t  REGION  | 

Fig.  5.7  Schematic  diagram  of  the  numerical  model 
(Burridge  and  Knopoff,  1967) 

Burridge  and  Knopoff  (1967)  also  presented  results  for 
a  one  dimensional  model  by  numerical  computation.  The 
computational  model  consists  of  10  masses  connected  by 


. 
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springs  and  driven  by  springs  coupled  to  a  moving  plate.  The 
1st  three  form  an  auxiliary  seismic  zone,  the  next  two  a 
viscous  zone  and  the  last  five  forms  primary  seismic  zone 
(fig.  5. 7). The  viscosity  is  so  high  for  the  middle  zone  that 
the  friction  force  is  never  exceeded  by  the  driving  force. 
Furthermore  dynamic  friction  is  assumed  to  be  a  function  of 
the  velocity  of  block  relative  to  the  friction  surface.  The 
assumed  friction  law  is 

ldKH 

_  6  —  -Ell  M.  ?  H 

I  +•  A  (u-h) 

B  _____  -  B  iX  ll 

=  1  -  ACci  +H) 

where  H  is  a  value  of  velocity  below  which  the  friction 
force  is  linearly  related  to  the  velocity.  E,  B,  A  are 
constants  incorporating  seismic  radiation  and  viscous 
effects . 

The  computation  shows  the  movement  of  the  6th  particle 
starts  first  followed  by  motion  of  the  7th  and  so  on  down  to 
the  10th  particle.  Then  motion  is  reflected  back,  particle 
10  imparts  further  motion  to  particle  9  and  so  on  in  the 
reverse  order  up  to  the  viscous  region.  This  event  may  be 
regarded  as  a  main  shock.  After  that  there  is  a  quiet 
period,  only  particles  4  and  5  move  according  to  a  slow 
adjustment  of  viscous  material  following  the  quake.  After 
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that 

1-3. 


the  mot 
This  quake 


on  of  particle  4  triggers  a  shock 
corresponds  to  an  after  shock  (fig 


in  region 
5.8)  . 


OUt  1  FI  RXX)  OUAKC 


VISCOUS  REGION 


J _ 1 _ 1 _ 1 _ L  -i  J  I  I  I  1  1  I  I  I  I  I  „ 

*9  \  2 *  *9  *9  O  C  J  2  5  4  .5  £ 

°  r>  «  r>  T»Wf 

o 


Pig  5.8 


For  a  number  of  shallow  focus  earthquakes  which  have 
surface  trace  of  rupture  it  has  been  found  that  there  exists 
a  correlation  between  the  magnitude  of  the  shock  and  the 


Fig.  5.9  Earthquake  magnitude  versus  fault  length. 

fault  length  (fig.  5.9).  King  and  Knopoff  (1968)  examined 
this  fact  in  a  laboratory  model  similar  to  the  Burridge  and 


Knopoff  (1967).  Assuming  that  the  drop  in  potential  energy 
of  the  system  is  a  measure  of  the  magnitude  of  the 
earthquake  and  the  number  of  blocks  moved  as  a  measure  of 
fault  length  (L),  they  plotted  the  data  and  found  that  there 
is  no  linear  relationship  as  observed  in  natural  events 
(fig.  5.10). 


Fig.  5.10  Shock  energy  versus  fault  length,  solid  dot  L=8; 
open  circle  L<8  (King  and  Knopoff,  1968). 


log  (L  d) 


Fig.  5.11  Shock  energy  versus  Ld  (King  and  Knopoff,  1968). 
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However,  better  linear  correlation  are  discovered  when 
UgE-  is  plotted  against  U^.LI>  ,  Lq  LT)^  ,  and  U^uX  .where 
X)  and  are  peak  and  average  displacement  respectively.  An 
example  of  the  correlation  between  B  and  LX  is  shown  fig. 
5.11. 


OtsuKa  (1972)  extended  the  Burridge  and  Knopoff 
computational  model  to  2D.  The  blocks  are  suspended  from  an 
overhead  support  by  leaf  springs  and  are  in  contact  with  a 
moving  floor.  There  are  2000  blocks  in  a  100X20  grid.  Here 
also  the  frequency  of  occurrence  of  event  and  energy 
released  (measured  by  number  of  blocks  moved)  shows  a 
relation  similar  to  the  Gutenberg- R i chter  relation  (fig 
5.12). 


Further  large  events  are 
regions  where  seismic  activity 
consistent  with  the  fact  that  in 


Fig.  5.12 

Magnitude-frequency  relationship 
N  is  the  number  of  events;  the 
horizontal  axis  gives  numbers  of 
blocks  that  slips  in  a  given 
cluster  (Otsuka,  1972) 

found  to  occur  in  those 
is  anomalously  low.  This  is 
the  period  preceeding  the 
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occurrence  of  a 
usual ly  low  in 
seismic  energy  is 


large  earthquake,  seismic  activity  is 
the  region  from  where  a  large  amount  of 
to  be  released  (fig.  5.13). 


During  mainshock 


Fig.  5.13  Darkened  area  shows  the  region  that  sliped 


(Otsuka,  1972). 


Dieterich  (1972b)  also  computed  earthquake  sequences 
with  a  one  dimensional  model  similar  to  Burridge  and 
Knopoff.  Four  cases  are  considered  in  which  the  coupling 


spring  and  friction 

are 

as  fol lows  ( 1 ) 

elastic 

spr i ng 

and 

time  independent  friction 

( 2 )  elastic 

spr i ng 

and 

t  i  me 

dependent  friction 

(3) 

vi scoe 1 ast i c 

spr i ng 

and 

t  i  me 

independent  friction 

(4) 

vi scoe 1 ast ic 

spr i ng 

and 

t  i  me 

independent  friction. 

The 

vi scoe 1 ast i c 

spr i ng 

is  taken  to 

have  some  partial  stress  recovery  mechanism  in 

order 

to 

produce  after-shocks. 

,  The 

time  dependent 

static 

friction  is 

taken  on  the  basis  of 

result  of  laboratory 

experiments 

and 

is  given  by  the  equation  5.1. 


89 


In  all  cases,  the  frequency  with  which  a  block  moves  is 
related  to  its  static  friction;  blocks  with  low  friction 
moves  more  frequently  than  those  with  high  f r i ct i on ( f ig . 5 . 1 4 


jS  OF  ELEMENTS 


Fig.  5.14 

Correlation  between  the  frictional 
strength  of  the  block  and  the 
number  of  times  the  block  slipped 
( D i e  ter i ch ,  1972b). 


) .  In  the 
even  when 


first  three  cases  no  after-shocks  are  produced, 
extreme  values  for  the  parameter  are  chosen  that 
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Plot  ot  the  earthquakes  in  the  model 
with  perfect  elasticity  and  simple 
friction.  There  is  no  distinction 
between  mainshock  and  aftershock 
(Dieterich,  1972b). 


give  large  stress  recovery  following  the  event.  For  example 
a  plot  of  the  first  case  is  reproduced  (fig  5.15).  There  is 
no  distinction  between  the  main  shocks  and  after-shocks.  The 
last  case,  viscoelastic  element  with  time  dependent 
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friction,  shows  an  after-shock  sequence  following  the 
main-shock.  These  occurs  on  the  same  segment  affected  by 
main-shock.  Further  the  duration  of  the  after-shock  directly 
related  to  the  magnitude  of  the  mainshock.  Large  earthquakes 
have  after-shocks  sequences  that  last  longer  than  those  of 
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Fig.  5.16 

Plot  of  earthquake  for  time 
dependent  friction  and  visco¬ 
elasticity.  Arrows  are  marked 
against  mainshocks  (Dieterich 
1972b  )  . 


smaller  earthquakes.  Also  the  frequencies  of  main-shocks  and 
after-shocks  are  easily  distinguishable  (fig.  5.16). 
Dieterich  (1973)  has  also  developed  numerical  schemes  for 
two-  and  three-dimensional  simulation. 


Cohen  (1977)  stimulated  earthquake  sequences  on  the 
basis  of  Dietrich's  (1972b)  model  and  concluded  that  seismic 
gaps  are  more  likely  to  occur  when  there  are  large 
variations  of  elastic  or  frictional  parameters.  If  the 
parameters  are  homogeneous  from  block  to  block  the  simulated 
earthquake  are  periodic  and  successive  events  propagate 
throughout  the  fault.  With  the  increase  in  the  spring 
constant  of  the  driving  spring  compared  to  the  connecting 
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spring,  the  length  of  shocks  decreases.  Furthermore  there  is 
a  correlation  between  the  displacement  of  the  largest 
after-shock  and  those  of  the  main-shock  (fig.  5.17). 
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Fig.  5.17  Plot  between  displacement  of  the  mainshock 
and  the  largest  aftershock  (Cohen,  1977). 

Recently  Rundle  &  Jackson  (1977)  have  re-examined  the 
Burridge  and  Knopoff  model  employing  20  blocks  with  elastic 
as  well  as  anelastic  coupling.  For  the  elastic  model  he  took 
three  different  friction  laws.  In  the  first  case  the  static 
friction  varies  violently  from  block  to  block  (maximum  to 
minimum  is  30:1)  while  dynamic  friction  lies  at  constant 
value  much  below  the  lowest  value  of  static  friction.  In  the 
second  case  static  friction  changes  very  smoothly  and 
dynamic  friction  is  0.8  times  static  friction.  In  the  third 
case  the  dynamic  friction  is  held  constant  and  static 
friction  changes  by  constant  amount  over  the  dynamic 
friction  from  block  to  block.  For  the  anelastic  model, 
static  friction  is  assumed  to  be  stress  as  well  as  time 
dependent  and  static  friction  -f'  decreases  at  a  rate 
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where  is  a  constant  and  s"  (t)  is  the  stress  acting  on 

the  block. 


In  all  these  cases  they  found  that  the 
frequency-magnitude  relationship  for  stimulated  earthquake 


MAGNITUDE  MAGNITUDE 


Fig.  5.18 

Frequency-magnitude  relationshi 
(Rundel  and  Jackson,  1977) 

( a  )  Elastic  mode  1 
(b)  Anelastic  model 


is  not  consistent  with  that  of  natural  earthquake  (fig. 
5.18).  With  the  different  value  of  elastic  parameter  the 
situation  does  not  improve.  So  they  have  concluded  that 
linear  behaviour  of  magni tude- frequency  relationship  is  not 
an  immutable  law  but  rather  is  dependent  on  the  mechanical 
properties  of  the  fault,  although  in  some  real  earthquakes, 
magnitude  frequency  relationship  is  similar  to  their  result. 
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King  (1975)  has  pointed  out  that  if  the  breaking  strength  is 
considered  to  be  finite  so  the  model  can  store  only  a  finite 
amount  of  strain  energy,  then  there  wi 1 1  be  a  decrease  in 
the  relative  number  of  large  event,  hence  the 
frequency-magnitude  relationship  may  not  be  linear. 

With  the  anelastic  model  it  is  possible  to  simulate 
aftershock  sequences.  Here  if  the  dynamic  friction  is 

l5cL 

greater  than  j  ,  the  stress  drop  during  the  first  rupture 
is  small,  and  the  stress  cr(t)  remains  large,  causing 
decrease  in  frictional  strength.  Since  friction  is  then 
reduced  while  stress  remains  high,  aftershocks  may  occur. 
Moreover,  the  aftershocks  will  not  be  confined  only  to  the 
region  of  primary  shock,  since  the  stress  induced  failure 
may  also  occur  in  a  region  adjacent  to  that  which  ruptured 
during  the  primary  event  (this  border  region  having  been 
stressed  by  the  elastic  elements  coupling  it  to  the 
di spl aced  region ) . 

On  examination  of  the  variation  in  the  frequency  of  the 
time  intervals  between  events,  some  nonPoisson  components  in 
this  distribution  of  interevent  periods  are  found.  Rundle 
and  Jackson  attributed  these  to  the  interaction  between 
adjacent  events.  The  motion  of  one  block  alters  the 
compression  of  the  spring  connecting  it  to  the  neighbouring 
block,  thus  altering  the  stress  and  the  time  at  which  the 
stress  overcomes  the  frictional  resistance. 
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None  of  the  models  discussed  so  far  have  accounted  for 
healing  or  stopping  of  the  slippage.  Moreover ,  due  to  the 
discrete  nature  of  the  models  there  is  fictitious  spatial 
variation  of  stress,  so  that  any  calculation  for  stopping  is 
unreliable  with  these  models.  A  realistic  picture  can  be  had 
if  a  continuous  model  is  considered.  Knopoff  et  al.  (1973) 
have  developed  a  continuum  model  by  letting  the  interblock 
distance  of  the  Burridge  and  Knopoff  discrete  model  to 
approach  zero.  They  thus  formulated  an  one-dimensional 
unilateral  rupture  problem.  They  showed  that  if  the  static 
friction  increases  or  the  difference  between  the  driving 
stress  and  dynamic  friction  decreases  sufficiently  in  the 
region  ahead  of  the  rupture  tip,  the  rupture  can  stop.  Thus 
some  variation  in  either  friction  or  stress  appears  to  be 
necessary  to  stop  the  rupture.  Israel  and  Nur  (1979)  have 
studied  the  effect  of  strongly  variable  stress  and 
frictional  strength  and  concluded  that  such  a  heterogeneity 
is  closely  linked  to  the  stopping  phase  of  the  fault  motion. 

Mikumo  and  Miyatake  (1978,79)  have  investigated  the 
earthquake  sequence  on  a  three  dimensional  block  model  with 
non-uniform  friction  and  spring  constants,  subjected  to  a 
time  dependent  shear  stress.  They  found  that  after-shock 
with  high  stress  drop  occur  in  and  around  the  unruptured 
part  of  regions  on  which  main-shock  has  already  occured, 
while  those  with  low  stress  drop  take  place  successively  in 
spaced  regions  so  as  to  fill  the  gaps  which  were  not 
ruptured  during  main-shock.  The  main-shocks  take  place 
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successively  in  adjacent  unruptured  regions  and  sometimes 
show  slow- speed  migrations.  Also  they  found  a  good  linear 
relations  in  logarithmic  scale  for  source  area  versus 
frequency  and  seismic  moment  versus  frequency  of  the 
generated  after -shocks . 


The  major  contribution  of  the  simulation  studies  has 
been  in  the  understanding  of  the  rupture  propagation  along  a 
fault,  in  developing  correlations  among  source  parameters, 
and  in  explaining  the  interplay  among  foreshocks, 
mainshocks,  aftershocks  and  creep.  It  appears  that  the 
simulation  studies  should  be  extended  to  two-  and 
three-dimensional  viscoelastic  systems,  possibly  with 
multiple  faults.  With  complicated  rheological  models,  the 
greater  number  of  adjustable  parameters  may  give  freedom  to 
simulate  completely  on  computer  a  natural  earthquake 
sequence.  This  may  be  useful  for  earthquake  predictions. 
Moreover,  there  is  no  direct  link  between  movements  of  the 
blocks  and  far  field  radiation  in  the  simulation  models 
developed  so  far.  So  a  synthesis  of  techniques,  such  as  a 
combination  of  a  sliding  block  description  of  rupture  with  a 
dislocation  theory  description  of  the  far  field  effect  might 
prove  fruitful . 


. 

■ 

6.  RECENT  TRENDS  IN  THE  SOURCE  THEORY 


In  the  recent  past  studies  in  the  source  theory  have 
been  directed  to  develop  more  realistic  models  by  taking 
into  account  the  layering  in  the  source  region,  effects  of 
the  free  surface,  and  anelasticity  of  the  medium.  Also  there 
have  been  attempts  to  account  for  the  energy  balance  during 
faulting. 

With  the  work  of  Lee  and  Teng  (1973)  on  the  radiation 
pattern  for  P  and  S V  waves  due  to  various  point  sources 
embedded  in  a  flat  multilayered  medium,  the  study  on  the 
effect  of  layering  started.  The  mathematical  tools  developed 
by  Ben-Menahem  and  Vered  (1973)  and  Helmberger  (1974)  where 
applied  by  Langston  and  Helmberger  (1975)  and  Bouchon  (1976) 
to  shallow  dislocation  sources  in  a  layered  half  space,  and 
they  obtained  analytical  expressions  for  the  displacements. 
To  synthesize  the  strong  ground  motions  produced  by  a 
propagating  dislocation  embedded  in  a  layered  medium,  Heaton 
and  Helmberger  (1977,1978)  modeled  the  source  as  a 
distributed  set  of  shear  dislocation  points  and  computed  the 
individual  responses  using  generalized  ray  theory. 
Bouchon ( 1 977 , 79 )  developed  an  exact  discretization  scheme 
for  the  elastic  fields  generated  by  such  a  fault,  and  the 
seismograms  obtained  were  then  combined  to  yield  the  space 
and  time  dependence  of  the  ground  motion  (Bouchon,  1980). 
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The  results  showed  the  strong  directivity  effects  on  the 
propagating  rupture.  Large  surface  displacements  and 
high-frequency  motions  were  confined  to  a  narrow  zone  around 
the  fault  and  to  the  region  which  extends  beyond  the  source 
along  the  trend  of  the  fault.  SH  and  Love  waves  were  the 
dominant  contributions  to  the  ground  shaking.  The  vertical 
displacement  were  small  but  exhibited  high  frequency 
content.  The  presence  of  low-velocity  surface  layer  had  a 
very  severe  effect  on  the  amplitude  and  duration  of  the 
ground  shaking. 

The  effect  of  a  free  surface  on  seismograms  has  been 
examined  by  Israel  and  Kovach  (1977)  and  Langston  (1978). 
They  found  that  depending  on  source  orientation,  depth,  and 
time  function,  errors  in  moment  and  corner  frequency 
determination  may  be  large  due  to  neglecting  the  free 
surface . 

A  few  efforts  have  been  made  to  study  earthquakes  in 
the  context  of  a  more  general  constitutive  relation  such  as 
Budiansky  and  Amazigo  (1976)  and  Burridge  (1977).  They  have 
assumed  a  model  in  which  the  vicinity  of  the  a  vertical 
fault  plane  is  represented  by  a  Maxwell  viscoelastic  half 
space . 

In  some  of  the  earthquakes,  permanent  elevation  changes 
in  the  vicinity  of  the  focus  have  occured.  This  suggests 
that  during  faulting  work  is  done  against  gravitational 
forces.  Several  authors  (Meister  et  a  1  .  (  1968)  ,  Jungels  and 
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Frazier  (1972))  have  noted  that  the  work  done  against 
gravity  in  a  dip-slip  faulting  appears  to  much  larger  than 
usual  estimates  of  energy  release  in  the  faulting 
( Gutenberg- R i chter  energy-magnitude  formula),  with  the 
implication  that  the  conventional  measure  of  the  earthquake 
energy  is  a  gross  under-estimate.  This  difficulty  was 
resolved  by  Dahlen  (1977),  who  defined  an  earthquake  by 
specifying  on  the  fault  surface  a  jump  discontinuity  in 
tangential  displacement  which  expands  within  a  uniformly 
rotating,  se 1 f -gr avi t a t i ng  earth  model. The  total  amount  of 
energy  released  by  such  an  idealized  earthquake  is  the  sum 
of  three  distinct  quantities:  kinetic  energy  of  rotation, 
gravitational  potential  energy  and  thermodynamic  elastic 
internal  energy.  Although,  the  total  energy  release  is  in 
general  considerably  smaller  than  any  of  its  three 
individual  constituents,  with  this  model  the  energy  release 
rate  is  simply 


where  is  the  ambient  stress  before  faulting,  <T"L'*  and 

j  J 

Ul  are  the  stress  change  and  displacement  associated 
with  the  faulting  and  integral  is  over  fault  surface. 

And  the  expression  for  the  work  done  against  gravity  is 

where  is  the  acceleration  due  to  gravity,  and  the 

integral  is  over  the  entire  volume  of  the  body. 
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Dahlen  pointed  out  that  although  the  gravitational 
energy  change  may  be  several  orders  of  magnitude  greater 
than  the  energy  release,  this  large  change  in  the 
gravitatioal  energy  is  approximately  compensated  for  by  a 
similar  change  in  elastic  energy  initially  stored  in  the 
Earth.  Thus  the  actual  energy  released  (i.e.,  seismic 
energy,  workdone  against  friction,  new  surface  energy)  is 
rather  small  difference  between  these  two  large  potential 
energy  changes . Savage  and  Walsh  (1978)  have  also  concluded 
the  same. 

The  present  state  of  art  in  the  source  theory  can  give 
us  the  seismic  motion  due  to  a  given  condition  of  stress  and 
material  properties,  but  in  order  to  have  a  realistic 
prediction  of  the  earthquake  occurrence  and  particularly 
foreshocks,  aftershocks,  and  other  related  phenomena,  the 
actual  distribution  in  space  of  tectonic  stress  field  and 
material  properties  in  the  epicentral  region  should  be 
considered.  And  here  is  the  challenge,  neither  those 
distributions  are  known  precisely  nor  analytical  techniques 
are  developed  to  tackle  such  a  heterogeneous  distribution. 
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